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I  INTRODUCTION 


I.  INTRODUCTION 


1. 1  Background 

Snow  accumulation  and  melt  are  an  Important  subsystem  within  the 
hydrological  cycle.  In  the  USA,  snowfall  and  snowcover  are  of  economic 
importance  (Lichfield,  1989),  with  possible  economic  hardship  resulting 
from  abnormally  low  seasonal  snowfalls.  Figures  1  and  2  show  that  in  the 
USA  the  mean  maximum  length  of  snowcover  duration  is  eight  and  a  half 
months  (September  15  -  June  1)  and  the  mean  minimum  duration  two  weeks 
(January  15  -  February  1),  with  only  seven  states  experiencing  less  than 
30%  of  the  sample  years  without  snowcover.  The  snowfall  area  extends  over 
many  different  climatic  zones,  each  with  a  characteristic  topography  and 
vegetation,  e.g.  tundra,  prairie,  high  mountain  and  desert.  These,  with 
their  accompanying  differences  in  snowcover  duration,  result  in 
characteristic  snow  environments.  A  detailed  knowledge  of  and  ability  to 
simulate  the  processes  occurring  in  these  snow  environments  are 
specifically  required  for: 

1)  Prediction  of  streamflows.  Predictions  are  usually  short-term  (1  day 
-  1  week),  concerned  with  flood  flows,  and  long-term  concerned  with 
seasonal  water  yield  for  domestic,  agricultural  and  energy  supply 
purposes . 

2)  Assessment  of  the  Impact  of  land  use  changes.  For  example,  these  can 
be  caused  by  silvicultural  changes  (Swanson,  1972,  Berris  and  Harr, 
1987,  and  De  Walle  et  al.  1977),  and  leisure  related  activities, 
Simons  (1988). 

3)  An  understanding  of  the  action  of  snow  and  snowmelt  as  a  pollutant 
carrier  and  concentrator,  e.g.  acid  snow  flushes,  Goodison  et  al . 
(1986). 

4)  A  research  tool  for  use  in  hydrology,  botany,  zoology,  ecology, 
geomorphology,  agricultural  studies  (e.g.  crop  protection)  and  as  a 
consideration  in  the  context  of  possible  global  climatic  change. 
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Figure  1:  Mean  date  of  snowcover  formation  (North  America) 
McKay  and  Gray  (1981). 


Figure  2:  Mean  date  of  snowcover  disappearance  (North 
America),  McKay  and  Gray  (1981). 
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This  report  is  concerned  with  the  development  of  a  model  to  reach  these 
requirements.  First  the  current  models  available  are  considered  in  the 
context  of  the  above  requirements  and  Che  objectives  and  scope  of  the 
model  to  be  developed  are  evolved.  The  snow  environment  C'l  be  modelled 
is  Chen  considered.  The  remainder  of  the  report  discusses  the  resultant 
model  and  its  development. 

1 . 2  Current  modellinB  methods 

This  report  is  primarily  concerned  with  Che  modelling  of  snowmelt  and 
therefore  this  review  will  consider  only  snowmelt  models,  or  those  which 
by  default  have  to  involve  snow  accumulation,  i.e.  those  in  which  the 
snow  system  is  modelled  throughout  the  season.  Unless  specified,  the 
term  'snowmelt  models'  will  Include  those  modelling  the  whole  season  and 
Involving  accumulation.  There  are  many  snowmelt  models  in  operation 
today,  lAHS  (1986).  Morris  (1985)  presents  a  classification  based  on 
the  operational  scale  (point  or  catchment)  and  the  melt  calculation 
method.  Classification  can  however,  sometimes  lead  to  confusion.  To 
avoid  this,  three  main  aspects  of  snowmelt  models  have  been  identified: 
snovmelt  model  development,  scale  and  melt  volume  calculation,  and  these 
are  discussed  below. 

1.2.1  Snowmelt  model  development 

The  wide  range  of  snow  environments  (table  1  and  section  1.1)  has  tended 
to  result  in  the  development  of  environment  specific  models,  especially 
chose  for  high  mountain  and  prairie  environments  rather  than  one 
all-encompassing  snowmelt  model  (Male  and  Gray,  1975).  The  historical 
development  of  snowmelt  models  mirrors  chat  of  rainfall-runoff  models, 
i.e.  developing  in  complexity  from  basic  equations,  e.g.  Darcy's  Law  to 
complex  distributed,  multiple  equations  models  (Anderson  and  Samblcs. 
1989).  Snowmelt  model  development  has  followed  three  paths: 

(1)  Models  developed  as  snow  accumulation  and  melt  subroutines  within 
larger  hydrological  models  which  call  subroutines  relating  to  all 
the  processes  of  the  hydrological  cycle,  e.g.  infiltration  and 
evapotransplratlon.  Pangburn  (1987),  for  example,  is  concerned 
with  the  development  of  the  snow  'box'  within  the  larger  MILKY 
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Table  1  :  Snow  environments  and  related  models 


Snow  environment 
type 

Pack  depth 

Vegetation  Melt  Example 

cover  characteristics 

Prairie/steppe 

Shallow 

Grassland 

cultivated 

Rapid 

Male  &  Gray 
(1975) 
Kuz'min 
(1961) 

High  mountain 

Deep 

plus  alpine 
permafrost 

Forest 

Prolonged 

Leaf  &  Brink 
(1973a  &  b) 

Temperate  lowlands 

Deep 

Mixed/ deciduous 
forest,  pasture, 
cultivated  mix 

Prolonged 

Dunne  &  Black 
(1971) 

Tundra 

Shallow  plus 
permafrost 

Tundra 

vegetation  - 
low,  sparse 

Rapid 

Everett  & 
Ostendorf 
(1988) 

i 
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Military  Hydrological)  model.  Morris  (1982)  considers  various 
models  for  Inclusion  into  the  SHE  (Systeme  Hydrologique  European) 
model . 

(11)  Independent  snowmelt  models,  e.g.  Anderson  (1976). 

(ill)  Detailed  process  specific  models.  These  are  concerned  with  aspects 
of  snowmelt,  e.g.  snow  metamorphism  or  capillary  water  movement 
(Colbeck,  1974),  These  models  are  physically-based  and  aim  to 
mathematically  simulate  the  environmental  process  as  accurately  as 
possible,  e.g.  Bohren  and  Barkstrom  (1979)  and  Colbeck  (1977),  The 
model  output  and  scale  are  specific  to  the  process  being  modelled 
(figure  3).  These  models  can  be  used  as  'building  blocks'  of  more 
complex  models  that  model  the  whole  melt  process. 

1.2.2  The  scale  of  snowmelt  models 

There  are  two  basic  scales  of  snowmelt  model,  those  that  simulate 
snowmelt  at  a  point  (Obled  and  Rosse,  1977)  and  those  that  simulate 
snowmelt  over  an  area,  l.e.  a  catchment  (Price  et  al.,  1976).  Point 
snowmelt  models  are  usually  extended  to  form  areal  snowmelt  models  by  a 
consideration  of  areal  snowfall,  this  and  other  variables  are  assumed 
homogeneous  over  the  area.  Ferguson  and  Morris  (1987)  consider  that  all 
areal  snowmelt  models  consist  of  a  basic  structure,  as  shown  in  tigure  4. 
There  are  four  sub-units  to  this  structure. 

(1)  A  meteorological  submodel.  This  enables  meteorological  data  from 
possibly  distant  point  sources  to  be  extrapolated  to  the  snowpack 
surface  (.as  in  the  case  of  remoter  catchment  sites).  Ferguson 
(1984)  uses  lapse  rate  to  calculate  air  temperatures  at  a  maximum 
height  of  1265  masl  from  a  meteorological  station  at  400  masl. 

Other  variables,  e.g.  precipitation  totals  and  evaporation,  can  be 
calculated  in  a  similar  way  (Reiners  et  al.  1984).  In  mountainous 
areas  temperature  and  precipitation  changes  with  altitude  will  be 
important.  In  other  environments,  e.g.  prairies,  other  factors, 
e.g.  wind  ( snow  drifting)  will  be  Important  spatial  variants.  A 
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Figure  3:  3e  Quervain  (1972;,  a  sohematic  diagram 
of  the  local  variation  of  temperature 
gradient  and  mass  flux  in  a  snow  struct 
Small  scale,  detailed  mass  flux  model. 


Figure  4:  Components  of  a  snowmelt  runoff  model,  Ferguson 
and  Morris  (1987) 


sophisticated  meteorological  submodel  would,  in  reality,  be  an 
areal  snow  accumulation  model.  In  many  areal  snowmelt  models 
initial  snowpack  extent  is  important.  Errors  in  the  estimation  of 
initial  snowpack  extent  can  manifest  themselves  in  the  depletion 
submodel  rendering  the  melt  output,  however  accurately  calculated 
by  the  snowmelt  submodel,  inaccurate. 

(11)  A  snowmelt  submodel.  This  calculates  the  volume  of  meltwater 

reaching  the  ground.  There  are  two  main  methods  used  to  calculate 
this,  the  index  method  and  the  energy-budget  method.  These  are 
discussed  further  below. 

(ill)  A  transformation  submodel.  This  routes  the  estimated  meltwater  to 
the  outflow  stream.  These  vary  in  complexity,  reflecting  the 
variation  in  current  hlllslope-runof f  modelling  techniques  (e.g. 
figure  5).  Some  transformation  submodels  are  Included  with  the 
snowmelt  submodel,  and  the  meltwater  movement  through  the  snowpack 
and  the  hlllslope  are  considered  together  (Dunne  et  al.  1976). 

(Iv)  A  depletion  submodel.  In  a  point  model  this  is  implicit  in  the 
snowmelt  submodel.  In  an  areal  snowmelt  model,  the  depletion 
submodel  is  important  as  it  effectively  attempts  to  describe  the 
pattern  of  snowmelt  and  the  resultant  snowcover  pattern.  Buttle 
and  McDonnell  (1987)  evaluate  five  such  models  (figures  6,  7  and 
8).  Depletion  models  are  in  an  early  stage  of  development,  and 
most  treat  the  areal  snowcover  they  represent  in  a  lumped  approach, 
1 .e .  the  percentage  of  bare  ground  may  be  calculated  but  there  is 
no  indication  of  the  distribution  of  that  bare  ground.  Unless  the 
depletion  model  is  operated  in  smaller  homogeneous  areas,  e.g. 
Buttle  and  McDonnell  (1987),  or  Leaf  and  Brink  (1977),  no 
indication  of  spatial  melt  patterns  can  be  obtaf.ed. 

1.2.3  The  calculation  of  melt-snowmelt  submodel 


Snownelt  models  can  be  classified  on  the  basis  of  the  method  they  use  to 
calculate  melt.  Older  models  used  simple  regression  techniques  (Hendrick 
and  DeAngells,  1976).  These  regressed  snowmelt  volume  against  a  number 
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H.O.F.  =  HORTON  OVERUND 
FLOW 


S.S.S.F.  =  SUBSURFACE 

STORMFLOW 


RETURN 

FLOW 


1  ‘ 

1  GROUND 

GROUND 

1  WATER 

WATER 

1 

1  RECHARGE 

FLOW 

I  j  =  INPUTS  ONTO 

SATURATED  AREA 


S.O.F.  =  SATURATION 

OVERLAND  FLOW 


=  SATURATED  AREA 


igure  5i  Meltwater  pathways,  Price  ar.Q  Hendrie 


OCPTH  Of  WATER  EOUIVALENT  M  SNOWPACR 


t.tl  •«.«  ••.t* 

%  OP  AREA  WTTH  W.E.$  OROMATE  VALUE 


Figure  6:  Schematic  representation  of  four  models  of  snowcover 
depletion:  (a)  spatially  uniform  water-equivalent, 

melt  at  snowpack  margins  (model  1);  (b)  spatially 

non-uniform  water-equivalent,  melt  at  snowpack  margins 
(model  2);  (c)  spatially  non-uniform  water-equivalent, 

spatially  uniform  melt  (model  3);  (d)  observed  pattern 

of  peaK  water-equivalent,  spatially  uniform  melt  (model  4). 
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of  variables,  e.g.  air  temperature,  wlndspeed  and  precipitation. 

However,  there  are  two  main  methods  of  calculating  snowmelt  volume  today: 

i)  Index  models.  In  these  models  meltrate  is  expressed  empirically  as 
a  function  of  one  or  more  variables.  The  most  commonly  used  index 
variable  is  air  temperature,  often  in  the  form  of  a  degree-day 
index.  Other  index  variables,  e.g.  vegetation  and  solar  radiation 
can  be  used  (Male  and  Gray,  1981  and  Martlnec  and  Range  1986)  can 
be  used.  An  example  of  a  typical  index  model  is  Ferguson  (1984). 
Melt  volume  is  calculated  by: 

Vi  =  Di  M  (1) 


where 

3  3 

melt  volume,  10  m 

degree  days  above  snowline  (degree-day  Index) 

M  constant,  nm  °C  ^  day  ^ 

Aj  snowpack  area ,  km^ 

where  is  given  by: 

DjS  m«^(0/y)-Nlwl(0/X)  HaoC^(0,;^)-MO/X^(0,fe)  (21 

♦(a-*) 

where 

X  minimum  air  temperature  in  the  24  hours  to  0900  GMT  on  day  1 

y  maximum  air  temperature  in  the  next  24  hours 

z  minimum  air  temperature  in  the  next  24  hours 

(see  figure  9) 

A^  is  calculated  using  the  areal  feedback  equation; 

Vi  *  -  ''i  w"'" 
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figure  9'-  Definition  diagram  for  calculation  of 
degree-days  above  0°C  at  snow  line. 


where 

2 

Initial  snowpack  area,  km 

initial  mean  water  equivalent  depth,  mm. 

It  can  be  seen  from  this  typical  example  that  index  models  do  not 
require  detailed  meteorological  inputs  and  they  can  operate  on  very 
basic  meteorological  and  snow  area  data.  This  is  an  advantage  in 
the  more  remote  areas  where  only  basic  meteorological  and  snow 
condition  coverage  exist.  Snow  area  data  can  be  collected  by 
remote  means,  i.e.  satellite  (Rango,  1988,  and  Merry  et  al.  1977), 
or  estimated  using  snow-area  elevation  curves.  Areal  index  models 
have  a  primary  output,  melt  volume  and  a  secondary  output,  snow 
covered  area.  Index  models  are  the  models  that  are  used  most 
commonly  for  operational  forecasting  purposes  (Roberge  et  al., 

1988)  and  are  also  most  often  used  as  the  'snow  input'  in  the 
larger  hydrological  models,  e.g.  Pangburn  (1987).  The  most 
advanced  index  model  is  Anderson's  (1973)  model.  This  bases  the 
calculation  of  melt  on  the  snowpack  energy  budget  but  uses  index 
methods  to  calculate  each  of  these  components.  For  example,  net 
radiation,  Q* ,  is  calculated  by: 


Q*  =  0.007  (T^  -  32) 


(4) 


where 

T^  ambient  air  temperature ,  °F 
Q*  net  radiation,  lnches/6  hr 

The  Internal  energy  exchanges  in  the  pack  are  not  modelled  on  a 
physically  based  criterion.  Index  models  do  not  attempt  to 
simulate  the  physical  processes  occurring  during  melt.  They  are 
concerned  with  the  manipulation  of  index  variables  to  produce  the 
correct  melt  volume. 
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(li)  Energy-budget  models.  Energy-budget  models  use  the  equation  for 
the  energy-budget  of  a  melting  snowpack  as  the  basis  of  the 
calculation  of  melt  and  they  use  physically-based  equations  to 
calculate  the  components  of  this  equation  (this  distinguishes  them 
from  the  Anderson  (1973)  type  of  index  model).  The  energy  budget 
equation  is  given  in  Male  and  Gray  (1981): 


<4 


0*  +  0  +0,  +  0  +  Q  -  dU/dt 
^  ^e  ^h  ^g  > 


(5) 


where 

0^  energy  flux  available  for  melt 

Q*  net  radiation  transfer 

latent  heat  (evaporation,  sublimation,  condensation)  transfer 
at  snow-air  interface 

Q,  convective  or  sensible  heat  transfer  at  snow-air  interface 
h 

Q  heat  introduced  to  pack  by  rain 

Qp  heat  Introduced  to  pack  by  rain 

dU/dt  rate  of  change  of  internal  (or  stored)energy  per  unit 
area  of  snowcover 


Meltrate  is  then  calculated  from 


Qm  = 


(6) 


where 


M 


p 


snowmelt  water  equivalent,  cm  day 
latent  heat  of  fusion,  kj  kg  * 

-3 

density  of  water,  kgm 

thermal  quality  or  the  friction  of  ice  in  a  unit  mass 
of  wet  snow 


Anderson  (1976)  has  a  very  strict  definition  of  energy-budget 
models  restricting  them  to  fully-distributed  physically-based 


16 


models  able  to  simulate  snovrpack  characteristics  (density,  thermal 
conductivity,  temperature  changes)  as  well  as  Incoming  radiation, 
etc.  Obled  and  Rosse  (1977)  and  Anderson  (1976)  are  examples  of 
these  true  energy-budget  models.  There  is  yet  no  catchment-scale 
fully  distributed  energy-budget  model.  Most  energy-budget  models 
are  point  models,  the  exception  being  Leaf  and  Brink  (1973).  This 
is  an  energy-budget  model  but  presents  a  very  simplistic 
calculation  of  the  snowpack  characteristics  (Anderson,  1973). 
Energy-budget  models,  because  they  are  physically  based,  are  much 
more  flexible  in  their  output,  i.e.  snowpack  temperatures, 
radiation  totals,  etc.,  can  be  obtained  as  output  in  addition  to 
meltwater  volume.  However,  they  do  require  more  input  data  in 
order  to  operate.  Anderson  (1976)  lists  a  minimum  of  five  Inputs. 

There  are  advantages  and  disadvantages  in  both  the  index  and  the 
energy-budget  methods  of  melt  volume  calculation.  Anderson  (1976)  tested 
his  energy  budget  model,  presented  in  that  report,  with  his  previous 
temperature  index  model  (Anderson,  1973),  figure  10.  He  concluded  the 
following : 

(i)  The  minimum  input  data  required  for  the  operation  of  an 

energy-budget  model  are  a  good  estimate  of  Incoming  solar  radiation 
and  measurements  of  air  temperature,  vapour,  pressure  and  wind 
speed . 

(ii)  Energy-budget  models  perform  better  than  index  models  when  applied 
to  open  areas,  because  of  the  greater  ability  of  the  energy-budget 
model  to  operate  under  variable  meteorological  conditions. 

(ill)  Under  the  stable  meteorological  conditions  of  the  forest  floor, 

index  models  give  similar  results  to  those  of  energy-budget  models. 

(iv)  Physiographic  factors,  i.e.  topography,  vegetation  cover  and 
climatic  factors,  should  be  considered  when  choosing  between 
energy-budget  and  index  models.  Generally,  the  more  varied  the 
physiography  and  climate,  the  better  the  performance  of  the 
energy-budget  model  when  compared  to  the  index  model. 
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covet  water  •  equivalent  (inni) 


(v)  When  extreme  conditions  need  to  be  modelled,  energy-budget  models 
should  be  used. 


Anderson  (1976)  also  considered  the  demands  upon  computing  time  of 
the  model  and  concluded  (in  1976)  that  the  energy  budget  model  was 
too  large  and  time  consuming  to  be  of  much  practical  use,  but  that 
with  some  modification,  not  at  the  expense  of  accuracy,  it  could  be 
made  more  manageable. 


There  are  therefore  many  snowmelt  models  in  operation  today.  In  order  to 
fulfill  all  of  the  requirements  for  snow  modelling  outlined  in  the 
introduction  (l.e.  streamflow  prediction,  environmental  changes,  research 
tool,  etc.)  a  snowmelt  model  must  be  of  the  energy-budget  type,  and  at 
the  catchment  scale.  There  is  a  need  for  a  catchment  scale, 
energy-budget  based  model  that  is  able  to  simulate  the  effects  of  slope 
angle,  aspect  and  vegetation  cover  impact  on  snowmelt  volume,  snowpack 
characteristics  and  subsequent  areal  snow  patterns. 


1 . 3  Objectives  and  scope 

The  aim  of  this  project  is  to  develop  a  catchment-scale  distributed 

snowmelt  model  for  areas  experiencing  seasonally  frozen  ground. 

The  model  that  has  been  developed  is  SNOMO  (SNOw  MOdel)  and  has  been 

designed  to  meet  five  basic  requirements: 

1)  the  primary  output  from  SNOMO  Is  the  spatial  pattern  of  snowcover 
over  a  catchment  and  subsequent  volume  of  snowmelt 

2)  SNOMO  is  based  on  the  energy-budget  approach  to  snowmelt  modelling 
and  is  designed  so  as  not  to  require  calibration 

3)  SNOMO  can  account  for  the  effects  of  spatial  variance  in  aspect 
and  slope  angle  on  snowmelt 

4)  SNOMO  can  account  for  the  effects  of  spatial  variance  in 
vegetation  on  snowmelt 
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5)  SNOMO  Is  flexible  enough  to  allow  for  choice  of  input  and  output, 
and  manipulation  of  environmental  conditions.  In  accordance  with 
the  spatial  nature  of  SNOMO,  output  is  graphical  as  well  as 
numerical 

Requirements  (4)  and  (5)  will  be  met  in  the  period  May  1989  -  May  1990. 

The  requirements  are  achieved  by: 

1)  development  of  the  FORTRAN-77  code  at  Bristol  University,  U.K. ,  on 
a  SUN-360  workstation 

2)  validation  of  SNOMO  using  field  data  supplied  from  and  collected 
at  Sleepers  River  Research  Watershed,  Danville,  Vermont. 

3)  evaluation  of  SNOMO  against  Leaf  and  Brinks'  (1973a  and  b)  existing 
spatial  snow  model 

1.4  The  snow  environment 


The  snow  environment  is  a  complex  system  of  inter-related  factors.  The 
primary  outputs  on  the  catchment  scale  are  snowmelt  discharge  and  the 
distribution  of  resld'ial  snowcover.  These  melt  outputs,  the  energy 
budget  components  that  result  in  melt  and  the  Internal  properties  of  the 
snowpack,  are  considered  in  general  terms. 

1.4.1  Residual  snowcover  pattern 

Snowpack  area  and  depth  vary  due  to  the  accumulation  pattern  of  the 
snowpack  (Itself  dependent  upon  factors  such  as  drifting,  topography, 
etc.)  and  the  variation  of  the  energy  budget  components  (variation  caused 
by  differing  slope  angle,  aspect,  vegetation  cover  and  the  prevailing 
meteorological  conditions).  Rau  and  Herrmann  (1982)  present  a  simplified 
version  of  this  spatial  and  temporal  variability  (fig.  11).  This  is  a 
graphical  representation  of  the  snowcover  distribution  with  elevatiai' 


A 


20 


HS  Snowdepth 

HSvv  Water  equivalent  of  snowdepth 

Cl  Critical  Interval,  meltwater  production 

is  initiated 

TI  Transition  Interval,  runoff  formation  is 

■nitiated 


Figure  Schematic  model  of  an  Alpine  snowcover 

store  • Hau  and  Herrmann,  'dS2; 
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r>''ar  the  snow  season  and  considers  a  plane  model  slope,  the  position  and 
elevation  of  which  are  defined  by  the  coordinates  of  X,  Y  and  Z.  X  and  Z 
are  the  axes  of  rotation  for  variations  of  slope  angles  and  azimuths. 

The  evolution  in  time  of  a  snow  profile  at  a  given  point  (Xy/Yy)  on  the 
slope  is  drawn  along  the  T-axls  (points  1,2  and  3).  Snow  profile 
variation  with  elevation  is  demonstrated  for  any  time  along 
the  y  axis  (points  4,2  and  5).  Figure  11  implies  differential  melt 
upslope.  The  pattern  of  snowcover  retreat  is  dependent  upon  the  scale  of 
the  area  over  which  melt  is  occurring.  There  are  three  general  scales t 

2 

1)  Large  (^lOkm  ),  e.g.  Rango  (1988).  Snowmelt  patterns  relate  to 
regional  meteorological  patterns,  e.g.  localised  precipitation 
and  major  altitudinal  changes,  e.g.  mountain  "snow-lines". 

Snowcover  is  usually  categorized  as  lOOX  or  OZ  cover. 

11)  Medium  (1-lOkm^),  e.g.  Leaf  and  Brink  (1973a  and  b) .  Snowmelt 

patterns  relate  to  major  topographic  features,  e.g.  hills  and  river 
valleys  (photo  1),  valley  orientation  and  vegetation  changes,  e.g. 
coniferous  forest,  pasture  boundaries. 

2 

ill)  Small  (^Ikm  ),  e.g.  Dunne  and  Black  (1971).  Snowmelt  patterns 
relate  to  smaller  topographic  features,  e.g.  hillocks,  drumlins, 
individual  hillslopes  (aspect  and  angle) ,  field  boundaries 
(drifting  and  shading)  and  Individual  objects,  e.g.  telegraph 
poles,  trees.  Snowcover  at  this  scale  is  very  patchy  and 
discontinuous . 

Therefore,  in  the  development  of  a  spatial  snowmelt  model,  the 
resolution  that  is  required  for  the  modelling  of  the  residual  snowcover 
is  very  Important. 

1.4.2  Snowmelt  discharge 

This  follows  a  diurnal  pattern  reflecting  the  decrease  or  cessation  of 
meltflows  during  the  night.  Snowmelt  is  facilitated  by  rain.  Rain 
during  melt  is  more  common  in  some  snow  environments,  e.g.  Oregon 
(Berrls  and  Harr,  1987)  than  in  others,  e.g.  prairies  (Male  and  Gray, 
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Photograph  1  :  The  drainage  network  is  highlighted  by  the 
distribution  of  the  snowcover,  which  is  concentrated  in  the 
channels,  Kameiin  and  Cook  (1967). 
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1975).  Flooding  during  snowmelt  Is  facillated  by  rainfall,  abnormally 
high  air  temperatures  and  the  topographic  and  vegetational  homogeneity  of 
the  catchment  (Hendrick  et  al.  1971).  Meltflow  occurs  through  the 
snowpack  and  on  the  surface  of  the  snowpack  and  as  throughflow  (through 
the  soil)  and  overland  flow  (over  the  soil  surface).  Snowmelt  can  cause 
soil  saturation  due  to  the  presence  of  an  Impermeable  frost  layer  in  the 
soil  at  depth  preventing  infiltration  or  rapid  melt  releasing  a  large 
volume  of  water  in  a  short  time  and  exceeding  the  Infiltration  capacity 
of  the  soil. 


1.4.3  Energy-budget  components 


The  energy-budget  of  a  melting  snowpack  in  the  absence  of  rain  is  shown 
in  figure  12  ,  Oke  (1987).  The  energy  balance  in  the  absence  of  rain, 
is  described  as: 


AQ,  "  Q*  +  Qh  +  Qe  '^g 

where 

A  Qg  net  heat  storage  term  (ie.  dU/dt) 

Q*  net  radiation 

sensible  heat  flux 
Qg  latent  heat  flux 

Qg  heat  Introduced  to  the  pack  from  the  ground 

A  Qjjj  latent  heat  storage  change  due  to  melting  or  freezing 


The  daily  energy  totals  relating  to  the  energy  balance  components  of 
figure  12  are  given  in  table  2.  Meltwater  runoff  rate  per  unit  area, 
^r,  is  calculated  by: 

.  lOOO  < 

l(.p 

whe  re 

Ar  meltwater  runoff  rate,  mm  day”^ 

Ji£  latent  heat  of  fusion,  Jkg”^ 

^  density  of  water,  kgra 
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ENERGY  FLUX  DENSITY  (W  m 


00  04  08  12  16  20  24 

TIME  (h) 


Figure  12:  Energy  budget  components  for  a  melting  snow  cover 
at  Bad  Lake  Saskatchewan  (51°N)  on  10  April  1978, 
Oke  (1987). 


Energy 


componencs 


(MJm  ^  day 


Q« 

2.02 

2.97 

-0.11 

<5h 

-0.84 

-^e 

-0.07 

0.11 

lilr  - 

8.92mmh 

Table  2  :  Daily  energy  corals  for  energy  componencs 
and  Che  derived  term  Ar  for  fig-  2 
(Oke,  1987) 


Table  3 

SELECTED  DAILY  E.NERCY  FLUX  TRAMSFER  (W/mV  DURING 
THE  MELT  PERIOD  IN  THE  ABSENCE  OF  VEGETATION  (BAD 
LAKE,  SASKATCHEWAN). 


Date 

(Dav/Mon/Yr) 

Q.n 

Qin 

Qn” 

Qh 

Qc 

Q« 

1  M/75 

8090 

-6320 

1770 

186 

-855 

-  45 

12-4/75 

9620 

-8480 

1140 

782 

26 

-  22 

14-4/75 

12290 

-9430 

2860 

13 

-395 

-  4 

17-3/76 

4630 

^500 

130 

1830 

-555 

64 

27-3/76 

7200 

-7720 

-  520 

1517 

-208 

-237 

2S-3/76 

7790 

-7120 

670 

70 

■201 

-II 1 

29-3/76 

9070 

-7660 

1410 

532 

-  60 

-180 

30-3/76 

9290 

-6040 

3250 

827 

140 

-270 

a 

poiiiive  values  indicate  an  ener|y  fam  by  the  snow. 

(he  daily  net  radiation  Hua  transfer:  ^  Otn  *  Qin- 

(Male  and  Gray,  1981) 
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Table  3  shows  some  comparable  results  for  other  years  at  the  same  station 
(Bad  Lake,  Saskatchewan).  A  discussion  of  the  energy- balance  of  a 
melting  snowpack  can  be  found  in  Oke  (1987).  The  energy-budget 
components  and  the  factors  that  affect  them  are  now  considered  in  turn. 

(1)  Net  radiation,  Q* 

Q*  is  composed  of  a  short-wave  component  and  a  long-wave  component 
and  can  be  described  as; 

Q*  =>  -  Kf  +  L|  -  Lt  (9) 

where 

incoming  solar  radiation 
K't  reflected  solar  radiation 
L  ^  incoming  longwave  radiation 
Lt  reflected  longwave  radiation 


(a)  Short-wave  component.  This  consists  of  incoming  (Kjr)  and 

outgoing  (Kt)  components.  K'J' can  be  divided  into  diffuse  and 
direct  components.  Diffuse  radiation  in  the  portion  of 
-  that  is  reflected  and  scattered  (i.e.  by  clouds)  reflected 
between  the  surface  and  the  atmosphere  (back-scattered),  and 
reflected  global  (direct  and  diffuse  KV)  radiation  from  the 
surrounding  terrain,  figure  13.  In  mountainous  terrain  the 
diffuse  radiation  input  is  affected  by: 

1)  decreased  sky  dome  due  to  surrounding  topography 
11)  angle  of  slope 

ill)  reflected  radiation  from  surrounding  surfaces 
iv)  anisotropic  distribution  of  diffuse  radiation  in  the  sky 
dome  (concentrated  nearer  the  solar  disc  and  the 
horizon) . 
v)  altitude 


4 


■A 
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Direct  beam  shortwave  radiation  Is  the  portion  of  K^that 
reaches  the  earth's  surface  without  being  absorbed  or 
diffused.  Direct  beam  radiation  on  a  slope  can  be  calculated 
using  spherical  geometry  if  the  slope  angle,  azimuth  (aspect), 
latitude,  date  and  hour  of  the  day  are  known.  Direct  beam 
radiation  also  increases  with  ai,.ltude  due  to  the  decrease  in 
optical  air  mass  with  the  decreased  thickness  of  the 
atmosphere.  Figure  14  shows  the  effect  of  slope  angle  and 
aspect  on  the  amount  of  direct  solar  radiation  received. 

Tables  4  and  5  and  figure  15  show  some  of  the  results  from 
Handler  and  Ishlkawa's  (1974)  calculations  on  the  effect  of 
slope,  aspect  ('exposure')  and  mountain  screening  on  the 
direct  beam  radiation  received  at  the  surface  of  the  McCall 
glacier,  Alaska  (69*^  19 'N).  Here,  screening  and  aspect 
effects  resulted  in  a  loss  in  duration  and  energy  content  of 
the  incoming  direct  solar  radiation.  On  a  slope  the  upper 
portion  will  receive  the  most  incoming  solar  radiation,  due  to 
its  longer  exposure  to  the  sun,  than  the  base  of  the  slope. 
This  is  demonstrated  in  Dunne  and  Black  (1971),  fig.  16), 

Snow  allows  some  transmission  of  incident  short-wave 
radiation.  The  decay  of  the  flux  into  the  snow  follows  an 
exponential  curve.  The  amount  of  short-wave  radiation 
reaching  any  depth  z  is  given  by  Beer's  Law: 

-CL?. 

K  (10) 


where 

u,  short  wave  radiation  incident  at  depth  z 
£,  base  of  natural  logarithms 
CL  extinction  coefficient,  m-1 

The  extinction  coef f icient ,  CL,  is  dependent  upon  the  nature  of 
the  transmitting  medium  (figure  17)  and  the  wavelength  of  the 
radiation.  The  depth  of  shortwave  penetration  can  be  as  much 
as  1  metre  in  snow,  Oke  (1987). 
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Jun«  21 


-2  -1 

Figure  Hi  Direct  solar  radiation  (cal, cm  hr  )  received 

on  different  slopes  during  clear  weather  at  50  N,  lat. 
Three  slopes  are  shown:  north,  south,  and  east-facing 
(west  would  be  a  mirror  image  of  east),  for  summer  and 
winter  solstice  and  equinox  (vernal  is  a  mirror  image 
of  autumnal).  The  lefthand  side  of  each  diagram  shows 
the  distribution  of  solar  energy  on  a  horizontal  surface 
(0  gradient).  The  righthand  side  of  each  diagram 
represents  a  vertical  wall  (90°  gradient).  The  top 
of  each  diagram  shows  sunrise  and  the  bottom  shows  the 
sunset,  Price  (1981). 


Tabu  4  Loss  in  duration  and  energy  for  the 
wHou  McCau.  Glacier  owing  to  the  screening 

EFFECT  OF  THE  SURROUNDING  MOUNTAINS  FOR  DIFFERENT 
SOLAR  DECUNATIONS 

Solar  declination  Loss  in  duration  Loss  in  energy 

deg  %  % 

23.5  32.6  10. 1 

20  39.9  13-4 

o  67.6  55.7 

—  10  93.6  87.8 


Table  5  Energy  received  on  the  McCall  Glacier 

FOR  different  SOLAR  DECUNATIONS 
A  horizontal  surface  without  any  screening  of  the  sun  is 
considered  to  receive  100%  energy 


Solar  declinalicn 
deg 

23-5 

Exposure 

Screening 

% 

Total 

Vr. 

/o 

98.7 

/o 

89-9 

to 

88.8 

20 

98-3 

86.6 

85.1 

0 

75-2 

44-3 

33-3 

— 10 

67.4 

12.2 

8.2 

31 


(a)  Midsummer,  solar ^ 
declination  23.5 


(b)  Equinox, solar^ 
declination  C 


I - 

1  49-49 
fTT\  49-49 

amm  49^99 


(c)  Autumn,  solar  ^ 
declination  -TO 


The  loss  in  energy  of  direct  solar  radiation  of 
McCall  Glacier,  Wendler  and  Ishikawa  (1974) 


Mareii  23,  1967  Morcn  24,  1967 


March  29,  1967  April  I,  1967 


Bore 


CZ]  Snow 


Figure  16:  Distribution  of  snow  cover  on  a  hillside  at 

noon  each  day  during  the  melt  period  of  March 
1967,  Dunne  and  Black  (1971). 
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Extinction  coefficient  (cm* 


Figure  17:  Extinction  coefl  ,ient  versus  snow  density,  Anderson  (1976) 


Figure  18:  Albedo  versus  snow  surface  density,  Anderson  (1976). 
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Reflected  shortwave  radiation  (Kf)  is  dependent  upon  the 
amount  of  incoming  shortwave  radiation  (KJ')  and  the  surface 
albedo  (Ot): 


Kf=  K4'(0£,) 


(11) 


Net  shortwave  radiation  (K*)  can  therefore  be  expressed  as: 

K*  =  Kid  -  DC)  (12) 

The  albedo  of  snow  varies  with: 

1)  Contamination  of  the  snow  by  dust  (leading  to  occurrences 
of  'red  snow'  in  the  Alps,  i.e.  contamination  by  red  Saharan 
dust),  forest  litter,  etc.  This  increases  with  the  age  of  Che 
snow  and  reduces  the  albedo, 

2)  Density.  This  is  again  related  to  the  age  of  the  snowpack 
and  the  amount  of  meCaraorphlsm  and  compression  it  has 
undergone.  Density  Increases  as  albedo  decreases  (fig.  18). 
Albedo  is  also  related  to  snow  grain  size  (which  is  related  to 
density),  Bergen  (1975). 

3)  The  nature  of  the  surface.  The  albedo  of  a  surface  will 
change  dramatically  when  the  snowcover  is  removed  or 
conversely  once  snowfall  occurs.  Price  (1988)  mentions 
overnight  albedo  changes  from  0.40  Co  0.90  at  Perch  Lake, 
Ontario  (Canada),  due  to  snowfall  on  deciduous  forest.  Figure 
19  demonstrates  a  similar  situation,  i.e.  snowcover  on  open 
fields  has  a  high  albedo  whereas  the  albedo  of  bare-tree  cover 
with  undersCorey  snowcover  is  much  lower.  Figure  20  shows 

the  close  relationship  between  point  albedo  and  snowfall. 

Density  and  contamination  changes  Increase  with  age,  as 
mentioned  above,  and  cause  a  reduction  in  the  albedo  of  snow. 
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0  1  2  3  4  5 

distance  (nm) 

Figure  I9i  Albedo  variation  with  vegetation  cover, 
O'Neill  and  Gray  (1972). 


FEB.  MARCH  APRIL 

DATE  — 

Figure  20:  Albedo  variance  in  response  to  snowfall, 
O'Neill  and  Gray  (1972). 


January  1971, 
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O'Neill  and  Gray  (1972)  summarize  these  changes  (fig.  21). 
Average  values  for  the  albedo  of  snow  are  found  in  Male  and 
Gray  (1981),  Oke  (1987)  and  Ballck  et  al.  (1981). 

Increased  absorption  of  K^will  occur  around  dark  objects 
contained  within  the  snowpack,  i.e.  tree  trunks,  leaves, 
twigs,  fence  posts,  etc.  This  is  because  of  their  lower 
albedo,  greater  thermal  conductivity  and,  therefore,  greater 
ability  to  absorb  shortwave  radiation  than  snow.  Localised 
melting  will  therefore  occur  around  these  objects. 


( b)  Long-wave  component 


This  consists  of  incoming  (L^)  and  outgoing  (Lt)  components. 
Net  longwave  radiation  (L*)  is: 

L*  -  L-i-  Lf  (13) 

L^  is  dependent,  in  the  absence  of  cloud,  upon  the  bulk 
atmospheric  temperature  and  emissivlty  (which  depends  upon 
distributions  of  temperature,  water  vapour  and  carbon  dioxide) 
following  the  Stephan-Boltzmann  Law.  Neither  atmospheric 
temperature  or  emissivlty  fluctuate  rapidly  and,  therefore, 

LV  Is  almost  constant  throughout  the  day.  Lf  is  also 
dependent  upon  temperature  (surface  temperature)  and 
emissivlty  (surface  emissivlty).  If  the  surface  is  a  full 

radiator,  l.e.fe  -  1,  then: 

0 

LK  erf'’’ 

0 


where 

ff*  Stephan  -  Boltzmann  constant 

T  surface  temperature 

0 
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V  BAD  LAKE(1971) 


0  1 _ I _ L  -  -1 _ 1 _ L  ■  .1 _ I _ L _ J _ 1_ 

0  2  4  6  8  10  12  14  16  18  20 

AGE  OF  SNOW  SURFACE  (days) 


Figure  21:  Variation  in  albedo  with  time  over  a  melting  Prairie 

snoupack  (Bad  Lake,  1971)  and  a  deep  mountain  snowpack 
(U.S.  Corps  o'"  Engineers,  1956),  O'Neill  and  Cray  (1972) 


However,  in  reality  ^  is  rarely  I  and  equation  (14)  becomes 

where 

surface  emlsslvlty 
p  atmospheric  emlssivity 


The  values  of  and  are  greater  than  their  atmospheric 
equivalents  and  because  varies  considerably  during  the  day 
Lt  la  larger  and  more  variable  than  .  L'f  has  a  maximum 
value.  This  is  when  the  snowpack  is  melting,  i.e.  =  0°C 
and  1  and,  therefore,  L ^  equals  316  Wm  the  black  body 

emlttance  at  0°C.  During  forest  snowmelt  1,  may  become  large 
because  the  canopy  absorbs  a  large  proportion  of  the  Kj/ 
incident  on  it,  heats  up  and  radiates  in  the  longwave  (Price 
and  Petzold,  1984).  This  can  result  in  a  large  L*  which, 
although  K*  at  the  forest  floor  may  be  small,  will  result  in  a 
large  Q*. 

More  detailed  accounts  of  the  net  radiation  exchanges  can  be 
found  in  Obled  and  Harder  (1978),  Oke  (1987)  and  Male  and  Gray 
(1981).  Q*  can  and  has  been  used  as  an  index  for  melt, 
instead  of  air  temperature. 

Price  (1988)  relates  Q*  to  and  in  order  to  use  Q*  as  an 
index  for  melt.  Olyphant  (1986)  table  6,  demonstrates  the 
relative  Importance  of  the  components  of  Q*  during  the  melt  of 
a  mid-latitude  alpine  snowpack.  Ambach  (1974)  discusses  the 
effect  of  cloudiness  on  Q*. 

11)  Sensible  and  latent  heat  exchanges,  and 

Tables  2  and  3  indicate  that  the  sensible  and  latent  heat 
exchanges  are  usually  of  secondary  Importance  when 
compared  with  net  radiation  during  snowmelt.  There  are. 
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Table  6 

Simulated  radiation  totals  for  individual  partial  areas  of  snov<  cover  in 
upper  Green  Lakes  Valley" 


Snowfield* 

K‘ 

Si 

(Ml  m-‘  d-') 

Dl 

(MJ  m-‘  d-) 

K 

■  (MJ  m-‘  d-) 

L, 

(MJ  m-‘  d.-') 

A 

(MJ  m-‘  d-')  - 

Arikaree 

11 

14.90 

7.70 

1.63 

19.38 

5.44 

8.37-17.50 

4.94-8.67 

0.80-3.68 

13.23-21.22 

2.85-13.90 

Green  Lake  5 

8 

17.20 

7.66 

1.63 

19.17 

5.65 

15.90-17.79 

7.28-7.91 

1.38-2.01 

18.33-19.63 

5.02-  7.03 

Kiowa 

8 

14.78 

8.04 

1.42 

19.88 

4.77 

13.06-15.91 

7.S7-S.25 

1.00-1.76 

19.38-20.89 

3.52-  5.53 

Rock  Glacier 

3 

12.01 

8.08 

1.59 

19.46 

5.11 

11.05-13.22 

8.00-8.20 

1.42-1,80 

19.30-19.80 

4.64-  5.40 

Green  Lake  a 

4 

16.03 

7.S7 

J.26 

19.97 

4.77 

15.70-16.41 

7.24-8.67 

0.38-2.05 

17.80-22.23 

1.63-  7.49 

Unobsinjcted 
ridge  top 

15.57 

9.21 

23.02 

‘Upper  entrj'  represents  the  mean  and  lower  entries  represent  the  range  of  simulated  irradiance  totals. 
•Individual  snowfields  are  identified  in  Figure  1. 

•Number  of  sites  evaluated. 


( Olypha-i  c ,  1 986 ) 
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however,  certain  environments  In  which  these  exchanges 
Increase  In  Importance,  e.g.  Golding  (1978).  Table  7 
gives  some  of  the  results  of  Golding's  study  of  the 
Importance  of  evaporation  during  Chinooks  In  Alberta 
(Canada)  and  demonstrates  that  In  most  of  the  sites 
Investigated  snowpack  loss  due  to  evaporation  was  greater 
than  that  due  to  melt. 

ill)  Ground  heat  flux, 

This  Is  a  negligible  energy-budget  component  when 
compared  to  the  others  (tables  2  and  3).  more 

Important  In  shallow  snowpacks  than  deep  snowpacks. 
affects  the  state  of  the  soil,  l.e.  frozen,  partially 
frozen  or  unfrozen,  which  affects  the  state  of  the  base 
of  the  snowpack  and  the  Infiltration  of  meltwater  during 
melt.  The  distribution  of  frozen  soil  will  vary  areally. 
The  forest  litter  present  under  deciduous  forest  protects 
the  soli  from  frost  penetration  prior  to  the  major  winter 
snowfalls.  Therefore  frost  penetration  Is  usually  much 
"^less  under  deciduous  forest  than  coniferous  and 

non-forested  areas.  Frost  penetration  In  the  open  is 
dependent  upon  the  timing  of  the  first  major  winter  snows 
and  the  preceding  meteorological  conditions.  If  the 
snowfalls  are  late  and  preceding  temperatures  are  very 
low,  then  frost  penetration  will  be  great.  If  the  soil 
Is  frozen  when  snowfall  occurs,  the  snow  will  tend  to 
freeze,  especially  If  at  the  start  of  winter  the 
temperatures  oscillate  around  zero,  or  there  Is  an  early 
winter  thaw  to  the  soil  resulting  In  a  snowpack  with  a 
frozen  base.  Included  in  Q  may  also  be  the  heat  evolved 
from  the  decomposition  of  organic  matter  present  between 
the  base  of  the  snowpack  and  the  soil. 
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Table  7  :  Mean  dally  snowpack  evaporation  and  melt  during  Chinooks, 
January  -  March  1975  and  1976.  (Golding,  1978). 


Site 

Elevation 

1975 

Dally  snowpack 

Evap  (mm)  Melt  (mm) 

1976 

Dally  snowpack 

Evap  (mm)  Melt  (mm) 

1 

1740 

0.9 

0.9 

1.9 

1.3 

2 

1700 

1.0 

0.6 

1.7 

1.5 

3 

1890 

0.8 

0.2 

1.6 

0 

4 

2190 

- 

- 

4.0 

0 

5 

2360 

2.4 

0 

- 

- 

6 

1280 

1.0 

0.5 

2.4 

2.4 

7 

1520 

0.9 

0.5 

2.1 

2.0 

8 

1690 

1.0 

1.7 

2.0 

3.2 

9 

1430 

1.2 

0.1 

- 

- 

10 

1580 

1.1 

1.6 

1.9 

5.5 

11 

1520 

1.4 

0.4 

1.8 

4.0 

12 

1520 

1.1 

0.5 

0.8 

3.7 

13 

1600 

- 

- 

0.8 

0.6 

14 

1120 

1.4 

0.1 

1.0 

3.9 

15 

1410 

1.0 

0.4 

0.8 

1.5 

16 

1630 

0.8 

0.2 

1.0 

0.4 

17 

2020 

1.0 

0 

2.2 

0 

18 

2200 

2.3 

0 

3.0 

0 

19 

2320 

0.8 

0 

3.8 

0 
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iv)  Raln-on-snow. 

Rainfall  on  snow,  commonly  called  rain-on-snow,  can  be  a 
contributary  factor  in  snowmelt  flooding  by  increasing 
meltrate  and  adding  to  total  discharge.  The  energy 
transferred  by  rain  to  the  snow  Is  the  difference  between 
the  energy  content  of  the  rain  directly  above  the  snow 
surface  and  Its  energy  content  upon  reaching  thermal 
equilibrium  within  the  back.  Two  situations  occur: 


a)  rainfall  on  a  melting  isothermal  pack  where  the  rain 
does  not  freeze 

b)  rainfall  on  a  pack  with  a  temperature  below  0°C, 
where  the  water  freezes  and  releases  its  latent  heat  of 
fusion  and  the  pack  warms  up 


v)  Latent  and  sensible  heat  storage  fluxes,  and 


represents  the  convergence  or  divergence  of  sensible 
heat  fluxes  within  the  snowpack,  i.e.  the  changes  in 
internal  energy  of  the  pack.  is  the  latent  heat 

storage  change  due  to  melting  or  freezing.  In  order  for 
the  melt  to  occur,  the  pack  first  has  to  reach  0°C  and 
then  the  latent  heat  change  can  occur,  i.e.  melt.  This 
is  demonstrated  In  figure  12.  In  the  morning,  the  heat 
input  into  the  pack  is  almost  entirely  put  into  storage, 
This  warms  the  pack  to  0°C  from  its  overnight 
temperature  (below  0°C).  Only  when  the  pack  is 
homogeneous  at  0°C  can  the  heat  input  into  the  pack  be 
used  to  change  the  proportions  of  ice  and  water  In  the 
pack,  i.e.  melt.  Melt  peaks  in  the  afternoon  and 
declines  as  Che  pack  cools.  This  cooling  is  retarded  by 
Che  meltwater  freezing  within  the  pack,  so  releasing  its 
latent  heat  of  fusion. 
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1.4.4  Internal  characteristics  of  the  snowpack 


Snow,  once  fallen,  undergoes  metamorphosis.  This  involves  a  change 
in  the  shape  of  the  snow  crystals  and  an  accompanying  change  in  the 
physical  properties  of  the  snow,  e.g.  density,  hardness,  thermal 
conductivity  and  thermal  diffuslvlty.  At  its  extreme  metamorphism 
is  responsible  for  the  conversion  of  fresh  snow  to  glacial  ice. 
Therefore  a  snowpack  composed  from  freshly  fallen  snow  (e.g.  an 
early  winter  pack)  will  have  different  physical  properties  from  one 
consisting  of  old  snow  (e.g.  a  spring  pack)  or  one  consisting  of 
layerrs  of  snow  of  variable  ages. 

Density  changes  result  in  a  compression  of  the  pack  height. 

Density  within  a  pack  is  usually  greatest  at  the  base  and  least  in 
the  surface  layers,  i.e.  density  Increases  with  depth.  Snow 
density  Increases  with  the  age  of  the  snow.  Fresh  snow  will  have  a 

_3 

density  value  of  around  0.l30gm  whereas  old  snow  exposed  by  melt 

-3 

will  have  a  value  around  0.400gm  .  The  density  of  the  top  of  the 

pack  can  be  increased  by  the  action  of  wind.  Figure  22a  shows  the 
effect  of  wind  on  increasing  the  density  of  the  top  of  the 
snowpack,  A,  Figure  22b  shows  the  same  profile  a  day  later,  after 
a  snowfall.  Layer  A  has  effectively  moved  down  the  pack  and  the 
fresh,  less  dense  snow,  sits  on  top  of  the  wind  affected  dense 
layer.  Layer  B  has  been  compressed  and  density  increases  as  depth 
Increases  overall.  Figures  22a  and  22b  also  demonstrate  the 
layering  of  the  pack  that  results  from  density  changes.  Figure  23 
presents  a  simple  layered  snowpack  where  the  'new'  snow  (less  dense 
is  easily  distinguished  from  the  'old'  snow  (more  dense).  Hardness 
Is  directly  related  to  density,  e.g.  ice  is  harder  than  fresh  snow. 

Figure  24  relates  the  change  in  density  to  the  change  In  effective 
thermal  conductivity.  Snow  density  is  directly  proportional  to 
snow  effective  thermal  conductivity.  The  thermal  conductivity  of 
snow  is  very  low,  which  makes  snow  a  good  Insulator  (of  houses, 
soil,  crops,  instruments,  etc,).  Thermal  conductivity  Is  also 
dependent  upon  the  temperature  and  micro-structure  of  the  snow. 
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Figure  22:  Density. versus  depth  measured  at  a  Whistler  Mountain,  B.C. 
ski  run  (elev.  1800  m) ,  Perla  and  Glenne  (1981). 


Moy  22,  1976 
SnowDQCk  Surface 


Density  (  g.  cm"^  Hordness  (kg.  cm*^) 


Figure  23!  Differences  in  snow  hardness  and  density  distinguish  new 
snow  from  old  snow  layers  within  the  snowpack, 

V/oo  and  Marsh  (1977), 
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Figure  24:  Approximate  relation  between  effective  thermal 
conductivity  of  snow  and  ice  and  density, 
Langham  (1981), 


The  thermal  diffuslvlty  of  snow  Is  also  low  and  directly  related  to 
density.  The  variations  of  snow  density,  thermal  conductivity  and 
thermal  diffuslvlty  that  can  occur  In  the  snowpack  result  In 
variations  of  the  snowpack  temperature  profile.  When  melt  occurs 
the  pack  has  a  homogeneous  temperature  profile  at  0°C  and  Is  said 
to  be  'ripe'.  Usually  temperatures  warm  from  the  top  to  the  base 
of  the  pack.  Inversions,  when  the  centre  or  base  of  the  pack  Is 
colder  chan  Che  top,  do  sometimes  occur.  Snowpack  layering  can 
also  occur  due  Co  the  refreezlng  of  melt  water  within  the  pack  and 
the  formation  of  extremely  dense  layers  (by  compaction  under  a 
large  snowpack  or  wind  action  at  the  surface).  This  results  In  Ice 
layers  or  Ice  lenses  within  the  snowpack. 

Therefore,  snowpack  physical  properties,  e.g,  density,  hardness, 
crystal  shape,  thermal  conductivity  and  thermal  diffuslvlty  change 
with  snowpack,  age  and  depth  due  to  metamorphosis.  These  changes 
result  in  compression  of  the  snowpack,  layers  within  the  snowpack 
and  a  variable  temperature  profile.  Metamorphosis  begins  as  soon 
as  the  freshly  fallen  snow  reaches  the  ground  and  increases  with 
time  and  snow  overburden  pressure.  Fresh  snow  and  old  snow  have 
characteristic  values  for  these  physical  properties.  Ice  layers 
can  form  because  of  metamorphosis  and  also  because  of  the 
refreezlng  of  melt  water. 

1.4.5  Summary 

Figure  25  summarizes  the  energy  exchanges  and  environmental  factors 
operating  In  an  Idealised  snow  environment  during  the  start  of  the  melt 
season.  Precipitation  falls  as  snow  or  rain  with  associated  snow 
accumulation  or  melt  enhancement.  Snow  accumulation  Is  not  areally 
homogeneous,  drifting  causing  greater  snow  depths  next  to  obstacles  (the 
fence)  and  at  the  base  of  slopes.  Snow  accumulation  will  also  be  less 
under  Che  forest,  especially  the  conifers,  than  In  the  open.  Snow 
Interception  also  results  In  albedo  variation.  This  and  other  factors 
result  In  the  variation  in  the  energy-budget  between  the  forested  and  the 
open  area.  The  slope  aspect  and  angle,  cloudiness,  wind  speed  and 
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other  factors  affect  the  energy-budget  for  the  whole  area.  Snowpack 
internal  properties,  e.g.  density,  thermal  conductivity  and  thermal 
dlffusivity  vary  with  age  of  the  pack,  pack  history,  l.e,  wind  hardened 
layers  and  pack  depth.  The  pack  can  be  multi-layered  and  contain  ice 
layers  and  lenses.  Melt  occurs  as  snow  surface  and  snow  throughflow  and 
as  throughflow  within  the  soil. 

The  residual  melt  pattern  on  the  small-scale  is  shown  in  figure  26. 
Snowcover  remains  in  areas  where  there  was  a  deeper  pack  initially,  i.e, 
at  the  slope  base,  surface  hollows,  and  by  the  fence,  where  shading 
retards  melt,  i.e.  at  forest  margins  and  where  the  energy-budget  is  less, 
i.e.  under  forest  cover.  This  diagram  demonstrates  the  variation  and 
hierarchy  in  the  importance  of  various  factors  in  Influencing  residual 
melt  patterns.  As  explained,  melt  is  usually  initiated  at  the  top  of 
slopes.  However,  in  this  case,  the  upper  part  of  the  slope  is  forest 
covered.  Field  observations  (photos  2  and  3)  have  shown  that  in  this 
situation  vegetation  cover  takes  precedence  over  slope  in  determining  the 
residual  melt  pattern.  Under  homogeneous  vegetation  cover,  photo,  slope 
determines  melt  pattern.  Under  different  snow  environments,  e.g,  the 
prairies,  other  factors  (e.g.  accumulation  patterns)  may  take  precedence. 
Snowmelt  flow  is  as  the  start  of  melt  situation  (figure  25)  but  the 
presence  of  bare  ground  and  saturated  soils  results  in  overland  flow. 

The  snowmelt  environment  is  therefore  a  complex  interrelated  system 
operating  under  variable  ground,  topographic,  vegetation  and 
meteorological  conditions.  The  effects  of  the  variance  of  these  factors 
on  the  energy-budget  of  the  snowpack  is  known  and  when  melt  occurs  there 
is  a  corresponding  residual  snowpack  distribution. 
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II  RESEARCH  DESIGN 
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II.  RESEARCH  DESIGN 


The  design  of  SNOMO  Involves  three  stages: 

(1)  development  of  the  model  structure,  i.e.  logic  pathways  and  source 
code 

(2)  validation  of  the  model  with  data  from  the  field 

(3)  evaluation  of  the  model  by  comparison  with  a  relevant  existing 
model 

The  emphasis  of  this  project  is  on  the  development  of  SNOMO,  therefore 
roughly  80%  of  the  project  time  is  spent  on  stages  (1)  and  (2)  and  the 
remaining  20%  on  stage  (3). 

2. 1  Development  of  model  structure 

The  basic  structure  of  SNOMO  is  shown  in  figure  27.  The  catchment  area 
that  is  to  be  modelled  is  subdivided  into  cells.  The  subdivision 
criteria  are  slope,  aspect,  vegetation  cover  and  altitude.  Each  cell  is 
assumed  areally  homogeneous  in  respect  of  these  four  factors.  Obviously, 
there  are  intra-cell  variations  and  these  Increase  with  the  size  of  the 
cell,  but  it  is  left  to  the  discretion  of  the  user  as  to  the  resolution 
required.  The  energy-budget  of  the  snowpack  in  each  cell  is  calculated 
at  the  mid-point  of  the  cell  (i.e.  mean  altitude).  The  melt  volume  for 
each  cell  is  then  calculated  by  multiplying  the  calculated  point  melt 
depth  by  the  area  of  the  cell. 

SNOMO  is  written  in  FORTRAN-77  and  developed  on  a  SUN-360  at  Bristol 
University.  The  program  is  divided  into  a  main  'control'  section  which 
governs  snowpack  characteristics,  melt,  and  handles  input  and  output  and 
a  large  subroutine,  TSTM,  which  calculates  the  energy-budget  components 
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and  the  temperature  of  the  snow.  The  subroutine  TSTM  is  adapted  from  the 
Terrain  Surface  temperature  Model  (Ballck  et  al.,  1981a).  Mr  Randy 
Scoggins  (WES)  helped  In  the  conversion  of  TSTM  into  a  subroutine. 

Appendix  i  shows  the  main  logic  structure  of  the  current  SNOMO. 

Individual  subroutines  and  equations  used  are  considered  in  Section  III. 

The  data  required  for  the  present  operation  of  SNOMO  is  given  In  table  8. 
This  appears  extensive  when  compared  to  that  required  for  an  Index  model. 
However,  the  meteorological  inputs,  with  the  exception  of  cloud  cover 
type  and  amount,  are  those  that  are  usually  measured  at  meteorological 
stations  even  In  the  remoter  sites.  Cloud  cover  type  and  amount  could  be 
estimated  from  synoptic  conditions  if  required.  Snow  absorptivity  (1  -  OL 
),  emlsslvlty,  dlffuslvity,  conductivity  and  density  all  have  recognised 
values  corresponding  to  their  age,  contamination,  etc.  The  remaining 
inputs  can  be  obtained  from  a  good  topographic  map  and  a  vegetation  map. 

If  a  vegetation  map  is  not  available  and  a  ground  survey  impossible,  then 
aerial  or  satellite  information  could  be  used. 

Subroutine  TSTM,  calculates  incoming  (direct)  solar  radiation  for  any 
slope,  angle  or  aspect.  This  makes  SNOMO  different  from  most  other 
snowmelt  models  which  either  require  a  measured  solar  radiation  input,  or 
calculate  solar  radiation  using  empirical  relationships  and  parameters. 
There  is  also  no  pre-melt  areal  snow  input.  This  is  because,  at  present, 
SNOMO  is  designed  to  accumulate  the  snowpack  from  pre-snow  conditions  in 
the  autumn  or  early  winter.  Accumulation  is  coarsely  modelled  using  snow 
precipitation  and  densities.  It  is  hoped  to  develop  SNOMO  to  enable  the 
initiation  of  modelling  on  any  day,  but  in  order  to  do  this  snowpack 
depth  for  that  day  will  need  to  be  known, 

2.2  Model  validation 


Model  validation  is  necessary  to  ensure  that  a  model  is  operating 
correctly,  l.e.  the  physical  equations  used  and  the  processes  they 
represent  are,  as  far  as  is  possible,  representing  reality. 
Validation  is  therefore  not  a  proof  or  disproof  of  the  utility  and 


54 


Table  8  :  SNOMO  operational  data  requirements 


Instrument  shelter  height 

Latitude 

Altitude 

Vegetation  cover  type 


For  each  cell : 

Air  pressure 

Cloud  cover  type 

Cloud  cover  amount 

Air  temperature 

Mind  speed 

Relative  humidity 

Surface  emlssivlty 

Surface  absorptivity 

Surface  moisture  content 

Thermal  dlffuslvity  at  depth 

Heat  conductivity  at  depth 

Slope  aspect 

Slope  angle 

Julian  day 

Precipitation  amount 

Snow  density 
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accuracy  of  a  model  but  a  useful  tool  in  its  development.  SNOMO  is  at 
present  a  prototype,  the  development  of  which  will  be  aided  by  validation 
data.  The  main  aim  of  this  project  (section  1.3)  is  to  develop  a  new 
model,  SNOMO.  Once  fully  developed  using  the  validation  data,  it  can 
then  be  'proved'  or  'disproved'  by  application  to  many  different 
catchments.  The  development  alone  will  take  three  years. 

The  data  required  for  the  validation  of  SNOMO  is  shown  in  table  9.  This 
data  is  detailed  and  there  are  few  snow  meteorlogical  stations  that  have 
adequate  data  coverage  of  this  type.  Validation  coverage  of  this  type, 
over  a  catchment,  is  currently  unavailable.  The  W3  catchment,  part  of 
the  Sleepers  River  Research  Watershed,  Danville,  Vermont,  USA  (figure  28) 
will  be  used  for  the  validation  of  SNOMO.  Snow  accumulation  depths  are 
not  as  great  in  New  England  as  elsewhere,  but  the  area  experiences 
variable  meteorological  conditions  during  melt  and  has  a  variable 
vegetation  cover,  both  of  which  present  a  good  test  to  the  modelling 
ability  of  SNOMO.  Snow  is  guaranteed  from  mid-November  to  late  March. 

2.2.1  Catchment  description  -  W3,  Sleepers  River  Research  Watershed 
Danville,  Vermont 


Research  began  at  W3  in  1965  with  the  initiation  of  a  NOAA-NWS 
cooperative  snow  research  project,  the  aim  of  which  was  to  improve  the 
data  base  available  to  snow  researchers.  The  main  meteorlogical  station. 
Townline,  began  operation  in  1968.  Data  from  Townline  enabled  Eric 
Anderson  to  develop  his  point  energy-budget  snow  model  (Anderson,  1976). 
The  data  from  Townline  is  comprehensive  and  reliable.  Research  at  W3  has 
continued  since  the  1970s.  At  one  stage  nine  snow  course  sites,  four 
basic  meteorological  sites,  and  three  weirs  were  in  operation  over  W3, 
now  only  three  meteorological  and  two  weir  sites  remain  in  operation. 
Research  at  '..'5  It  now  conducted  by  CRREL  (US  Army  Cold  Regions  Research 
and  Engineering  Laboratory,  Hanover,  New  Hampshire)  Mr  Hugh  Greenan  and 
Dr  Timothy  Pangburn,  and  the  UVM  (University  of  Vermont),  Mr  Bill 
Roberts,  UVM  initiated  a  ground  temperature  project  in  1985.  A  series 
of  frost  tubes  measuring  the  depth  from  the  surface  of  frozen  ground,  is 
in  operation  under  various  vegetation  cover  types  over  W3.  Weekly 
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Table  9  :  SNOMO  validation  data  requirements 


Net  shortwave  radiation 
Net  longwave  radiation 
Sensible  heat  exchange 
Latent  heat  exchange 

Areal  distribution  and  depth  of  the  snowpack 

Volume  of  snowmelt  and  river  discharge  over  the  snowmelt  period 
Density  of  the  snowpack 
Surface  emisslvity  and  albedo 
Precipitation  type  and  amount 


Snow/soll  saturation,  thermal  diffusivlty  and  heat  conductivity 
Soil  profile  temperatures 


Any  additional  information,  e.g.  snowpack  stratification,  snow  profile 
temperatures,  or  unusual  meteorological  conditions. 


•  Hawaii 


Vermont 


SLEEPERS  RIVER 
RESEARCH 
WATERSHED. 
DANVILLE.  VERMONT 


J  y,  45- 

u.s.A. 


'  Sleepers  River  j 

I  Waters  tied  / 

J  •Burlington  y 

1  0  StJohnsbury 

Montpelier  j  I 


VERMONT 


/CRREL 

'(Hanover] 


SUBCATCHMENT 

W3 


Figure  28:  Location  of  the  W3  watershed. 


readings  are  made  of  the  frost  depths. 

The  W3  catchment  is  considered  representative  of  the  land  use  cover,  soil 
and  topographic  conditions  found  in  the  northern  areas  of  the  New  England 
States  and  southern  Quebec.  Photograph  Q,  is  typical  of  this  New  England 
environment,  rolling  wooded  hills,  with  small  areas  of  pasture  land.  A 
more  detailed  description  of  W3  is  given  in  Anderson,  et  al .  (1977)  and 
Pionke,  et  al.  (1978). 

2 

The  basin  has  an  area  of  8. A  km  and  varies  in  altitude  from  1135  to  2280 
feet  above  mean  sea  level  (figures  29  and  30).  Vegetation  cover  is 
predominantly  forest,  coniferous,  deciduous  and  mixed  with  some  areas  of 
open  pasture  (figure  31).  There  is  no  arable  land  at  W3.  The  main 
deciduous  species  are  Birch  (yellow  -  Betula  allegheniensis ,  white  -  B. 
papyrifera  and  grey  -  B.  populifolia.  Beech  (Fagus  americana)  and  Maple 
(sugar  -  Acer  saccharum  and  red  -  A.  rubrum).  The  major  coniferous 
species  are  Red  Spruce  (Picea  rubra)  and  Balsam  Fir  (Abies  balsamea) . 

There  has  been  and  continues  to  be  forestry  activity  in  certain  areas  of 
W3,  mainly  in  the  coniferous  areas,  which  has  resulted  in  large  tracts  of 
clearcut,  A  lot  of  this  activity  occurred  in  the  early  1980s,  with  the 
result  that  any  new  tree  growth  in  these  areas  is  still  relatively  young 
and  light-loving  shrubs,  e.g.  Dogwood,  tend  to  predominate. 

The  Walts  River  formation  forms  the  underlying  solid  geology  of  W3.  The 
Walts  River  formation  consists  of  calcareous  granulites  and  calcareous 
schists  (composed  of  quartz,  calcite,  muscovite  and  biotite)  ,  quartz  mica 
schist  and  minor  micaceous  quartzite.  Superficial  deposits  are  a  till 
with  a  clay-silt  matrix  with  small  patches  of  poorly-sorted  horizontally 
stratified  (lenticular  bedded)  and  crossbedded  sand,  pebbly  sand  and 
pebble  gravel. 

Soils  are  in  general  fairly  shallow  and  are  variable  in  depth  over  small 
areas,  i.e.  from  zero  to  several  Inches  over  exposed  rock  outcrops  to  A 
to  6  feet  between  these  outcrops.  The  distribution  of  soils  as  mapped 
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Figure  29:  Topographic  map  of  the  W3  watershed, 
Anderson  et  al.  (1977) 


ELEVATION-FT.(MSL) 


PERCENT  OF  AREA  ABOVE  GIVEN  ELEVATION 


Figure  30:  Area-elevation  curve  for  the  W3  watershed 
Anderson  't  al.  (1977). 
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and  classified  by  SCS  Soil  Surveys  1960-64,  compiled  in  1969,  (Pionke  et 
al , ,  1978)  and  the  dominant  slope  angles  of  W3  (from  the  same  source  as 
the  soils)  is  shown  in  figures  32  and  33.  Table  10  shows  the  percentage 
areas  of  Hydrologic  Soil  Groups,  Land  Use  and  Slope  Distribution  for  W3 
in  1978  (Pionke  et  al.,  1978).  It  is  probable  that  the  figures  for  the 
Hydrologic  Soils  Group  and  the  slope  distribution  have  changed  relatively 
little.  Land  use  has,  however,  changed,  there  now  being  no  arable  land 
and  a  decrease  in  forest  due  to  clearcutting. 

2. 3  Model  evaluation 

It  is  necessary  to  compare  the  results  from  SNOMO  with  the  results  from  a 
model  which  most  closely  resembles  SNOMO,  as  an  indication  of  the 
relative  performance  of  SNOMO.  SNOMO  will  be  evaluated  against  Leaf  and 
Brinks'  (1973a  and  b)  model,  WATBAL,  figure  34.  In  a  review  of  models 
developed  for  forest  hydrology  (USFS,  1980),  WATBAL  was  selected  as  the 
most  "readily  useful",  along  with  PROSPER  (a  non-snow  environment  model). 
The  alms  and  model  structure  of  WATBAL  most  closely  resemble  that  of 
SNOMO  because  WATBAL  is  specifically  a  catchment  model.  Anderson's 
(1976)  point  model  most  closely  resembles  the  radiation  computations  and 
layering  structures  used  in  SNOMO,  but  it  is  a  point  model  and  does  not 
allow  for  vegetation. 

WATBAL  models  (1)  winter  snow  accumulation,  (2)  the  energy  balance,  (3) 
snowpack  condition,  and  (4)  resultant  melt  in  time  and  space  under  a 
variety  of  conditions.  WATBAL  can  model  the  effects  of  three  watershed 
management  practices: 

1)  clearcutting  in  patches 
il)  selective  cutting  (thinning) 
ill)  weather  modification  (cloud  seeding) 

Figure  35  is  a  diagrammatic  comparison  of  SNOMO  and  WATBAL.  There  are 
similarities  between  the  two  models.  They  are  both  energy-budget  models 
and  subdivide  the  catchment  into  computational  cells.  Differences  in  the 
two  models  are  listed  in  table  11. 
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Source  ;  PIONKE  el.al.  (1978) 


Source  ;  PIONKE  el.al.  0978) 


Table  10:  Percentage  areas  of  Hydrologic  Soil  Grouos ,  band  Use 
and  Slope  Distribution  over  W3  (Pio’-ke  et  al .  ,  1978  ) 


Parameter 


Mapping  Unit 


“  Area 


Hydrologic  Soils  A 

Group*  3 

r 

D 


0 

3 


Land  Use 


Cultivated 

Forest 

Pasture 

Idle 


67 

19 

3 


Slone  Distribution 


0-3% 
3-8 
?-l  5 
15-25 


38 


’  ? 


★ 


A 


(low  runoff  potential).  High  infiltration  rates.  Well 
excessively  well  drained  sands  or  gravels.  High  race  o 
water  transmission. 


3  =  Moderate  infiltration  rates.  Moderate  race  of  water 

transmi ssion. 


C  =  Low  infiltration  races.  Slow  rate  of  water  transmissit 

0  =  (High  runoff  potential).  Very  slow  infiltration  rates. 

Clay  soils,  soils  with  a  permanent  hith  water  table,  so 
with  a  claypan  or  clay  layer  at  or  near  the  surface  and 
shallow  soils  over  nearly  impervious  .material.  Very  si 
rate  of  water  transmission. 
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Table  11 


Comparison  between  SNOMO  and  WATBAL 


Variable 


SNOMO 


WATBAL 


Application 


Catchments  with  a  wide 
range  of  vegetation 
cover,  elevation  and 
meteorological  cond¬ 
itions  . 


The  North-American  sub-alpine 
environment,  i.e.  high  elevations 
with  almost  continuous  coniferous 
tree  cover  are  relatively  stable 
meteorological  conditions.  An 
application  of  WATBAL  to  the  Wolf 
Creek  watershed  (Oregon)  showed 
problems  in  applying  WATBAL  to  a 
sub-alpine  area  with  a  maritime 
climate  (USFS,  1980). 


Catchment 

division 

Meltwater  is  routed 
through  soil  and  there¬ 
fore  the  position  of 
the  cells  is  uncon¬ 
strained.  Variable 
size  of  cell. 

Meltwater  is  not  routed  through 
the  soil,  therefore  the  cells  must 
border  on  a  stream.  Cell  size  is 
an  average  of  10%  of  the  catchment 
area. 

Calculation  of 

See  section  III.l 

Generated  empirically  as  a 

shortwave 

function  of  maximum  temperature 

radiation 

and  slope/aspect  characteristics. 

Calculation  of 

See  section  III.l 

Empirically  adjusted  for  the  snow 

longwave 

environment  with-' without  cloud 

radiation 

cover 

Calculation 
of  cell 
precipitation 
values 


Cell  precipitation 
values  are  assumed 
identical  to  that  of 
the  nearest  meteorol¬ 
ogical  station. 


Base  station  is  at  a  lower 
elevation  than  the  catchment. 
Base  station  precipitation  is 
adjusted  until  the  specified 
peak  water  equivalent  at  each 
station  Is  reached.  Requires  a 
knowledge  of  the  peak  WE  of  each 
cell  and  implies  retrospective 
modelling 


Calculation  of 
cell  temperature 
values 


Cell  temperatures  are 
assumed  identical  to- 
that  of  the  nearest 
meteorological  station 


Base  station  is  at  a  lower 
elevation  than^t^egCa^chment . 

subunit  base 
where  A  and  B  are  empirical 
coefficients. 


Snowpack  Maximum  of  6  layers,  The  snowpack  is  treated  as  a 

layering  each  layer  with  single  homogeneous  layer, 

individual  properties. 
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A  copy  of  MATBAL  and  an  operational  data  set  have  been  obtained  from 
Charles  Troendle,  U.S.  Forest  Service,  Fort  Collins,  Colorado. 
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Ill  MODEL  DETAILS 


III.  MODEL  DETAILS 


3. 1  Introduction 

SNOMO  consists  of  a  'control'  program  and  various  subroutines,  figure  I 
(appendix  D-  This  section  considers  the  various  elements  of  figure  i 
in  detail  as  follows; 

(a)  snowpack  characteristics  and  structure. 

(b)  data  input  and  initialization  of  the  model. 

(c)  calculation  of  the  snow/soll  energy  budget  variables,  snow 
surface,  and  internal  temperature 

(d)  raln-on-snow 

(e)  calculation  of  meltrate 

(f)  snowfall  and  compaction  changes  to  pack 
characteristics 

(g)  melt  changes  to  pack  characteristics 

(h)  output 

3.2  Snowpack  characteristics  and  structure 

The  snowpack  is  modelled  as  either  a  2-layered  or  1-layered  system.  The 
snow  is  described  as  either  'old'  or  'new'  (fresh  snow).  The  single 
layered  pack  is  therefore  either  all  old  snow  or  new  snow  on  top  of  old 
snow.  The  pack  characteristics  and  temperature  are  held  in  two  matrices, 
DEPMX  and  THKMX.  There  are  four  matrices  in  total,  DEPMXl  (DEPMXl  for  1 
layer  pack),  DEPMX2  (DEPMX  for  2  layer  pack),  THKMXl  (THKMX  for  1  layer 
pack)  and  THKMX2  (THKMX  for  2  layer  pack),  see  fig.  36.  The  matrices  are 
also  defined  by  two  variables,  NOMATL  (number  of  layers)  and  NIPTS 
(number  of  Iteration  points,  i  .e .  the  layer  boundaries),  as  detailed 
below: 


72 


73 


i)  DEPMX  (depth  matrix) 

DEPMXl  Is  a  2  X  2  matrix  and  DEPMX2  Is  a  3  x  2  matrix 
1st  column:  depth  of  the  layer  boundaries,  centimetres 
2nd  column:  corresponding  temperature, 


11)  THKMX  (thickness  matrix) 

THKMXl  Is  a  1  X  6  matrix  and  THKMX2  Is  a  2  x  6  matrix 
1st  column:  thickness  of  the  layer,  centimetres 


2nd  column: 
3rd  column: 
4ch  column: 
Sth  column: 
6ch  column: 


2  —  1 

thermal  dlffusivlcy,  cal  cm  min 

2  —1  “1 

heat  conductivity,  cal  cm  min  k 
surface  emlsslvlty,  decimal 
surface  albedo,  decimal 


density,  gm 


-3 


The  5  physical  variables  In  (11)  above  refer  to  the  snowpack  and 
vary  with  the  age  and  depth  of  the  snowpack.  Fresh  snow  and  old 
snow  have  characteristic  values  (cable  12).  AC  present  Che  extreme 
values  for  old  and  new  snow  (i.e.  old,  dirty  and  end-of-melc  snow 
and  freshly  fallen,  clean  snow)  are  used  In  SNOMO.  These 
characteristics  are  held  In  two  vectors,  VNEWSN  (fresh  snow  values) 
and  VOLDSN  (old  snow  values).  These  vectors  are  1x6  and  the 
first  colimtn  Is  0.0.  This  Is  where  Che  thickness  of  the  pack  is 
held  in  THKMXl  and  THKMX2.  The  remaining  5  columns  are  the  same  as 
THKMXl  and  THKMX2.  The  vectors  VOLDSN  and  VNEWSN  are  substituted 
Into  the  thickness  matrix  when  and  where  appropriate.  There  Is  a 
third  vector,  VSOIL,  which  contains  the  same  variables  as  VNEWSN 
and  VOLDSN  but  for  the  soil  (see  table  12)  VSOIL  Is  used  when  snow 
disappears  due  to  melt  or  compaetloo»  or  when  SNOMO  is  at  the  start 
of  a  run  and  Che  snow  has  yet  Co  fall.  The  values  used  In  VSOIL 
depend  upon  soil  type;  If  this  Is  known  then  the  values  are 
obtainable  from  various  sources  such  as  those  used  for  snow. 


There  Is  a  threshold  snowpack  depth,  Che  critical  depth  of  5cm 
which  demarks  Che  presence  an  !  absence  of  the  snowpack.  A  depth  of 
5cm  was  chosen  as  snowpack  uepChs  under  5cm  arc  difficult  to 
measure  accurately  and  Che  loss  of  less  Chan  5cm  of  snow  to  the 
water  Input  of  SNOMO  was  deemed  negligible,  or,  In  the  case  of 
melt,  could  always  be  added. 
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Table  12  ;  Values  for  the  physical  characteristics  of  snow 


New  snow 

Old  snow 

Sandy  soil 

Clay  soil 

Density 
(gm  ) 

0. 10 

0.48 

1.60'^ 

1.60'*' 

The rmal 
dlf ^uslvljy 
(cm  min  ) 

0.06 

0.24 

0.36* 

0.39* 

Heat  0.03 
conduct^vity_j^ 

(calcm  min  K  ) 

0.08 

0.08 

0.09* 

0.13* 

Emlsslvity 

(decimal) 

0.82 

0.99 

0.91-0.93* 

0.88-0.97* 

albedo 

(decimal) 

0.95 

0.40 

0.43-0.33* 

0.60* 

Sources  : 


Oke  (1987),  Ballck  et  al . 
and  Gray  and  Male  (1981) 


( 1981a) 


3.3  Data  Input  and  model  Initialization 


Data  is  input  into  SNOMO  in  a  file  SNOMO.DAT  and  is  also  present  in  DATA 

statements  in  the  SNOMO  code.  The  file  SNOMO.DAT  contains  physical  and 

computational  data  and  variable  daily  data,  see  table  13.  The  data  file 

SNOMO.DAT  is  considered  below: 

LINES  1-3  :  discussed  in  section  3.2 

LINE  4 

Cl  Air  pressure;  This  is  input  in  millibars  and  measured  on  site 
or  taken  from  synoptic  charts 

C2  Cloud  type:  Index  values  ranging  from  1  to  8  are  used  to  determine 
cloud  type,  see  table  14.  If  cloud  type  data  is  unavailable,  then 
it  is  estimated  from  known  synoptic  conditions,  precipitation,  etc. 

C3  Observation  instrument  height:  This  is  the  height  of  the 

meteorological  Instruments  from  the  ground  or,  in  Che  case  of  wind 
speed,  from  Che  snow  surface. 

LINE  5 

Cl  Slope  angle:  Angle  of  the  slope  from  the  horizontal  measured  in 

degrees.  This  is  assumed  homogeneous  for  each  cell  and  can  either 
be  measured  in  the  field  or  estimated  from  topographic  maps. 

C2  Slope  aspect:  Measured  in  degrees  from  south  with  positive  values 
increasing  westwards  to  north  and  negative  values  increasing 
eastwards  to  north,  i  .e .  south  =  0°,  north  =>  I80'’,  due  east  ■  -90° 
and  due  west  »  +  90°.  Slope  aspect  is  assumed  homogeneous  for  each 
cell  and  can  be  measured  in  the  field  or  estimated  from  topographic 
maps . 

C3  Start  date:  The  date,  Julian  calendar,  of  the  first  day  modelled. 
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Table  13  ;  SMOMO.DAT. 


LINE  1 ,  VNEWSN 

LINE  2,  VOLDSN 

LINE  3,  VSOIL 

Cl  0.0  , 

C2  therruTl  diffusivity  (ctn  ^  ^ 

C3  heat  conductivity  (cal  cm  “  min  K  M 

C4  surface  emissivity  (decimal) 

C5  surface  albedo  (decimal) 

C6  snow  density  ( VNEWSN, _VOLDSN,  gm  ■^) 

soil  density  (VSOILgm  ) 

LINE  4 

Cl  air  pressure  (mb) 

C2  cloud  type  (index) 

C3  observation  instrument  height  (cm) 

LINE  5 
Cl 
C2 
C3- 
C4 

LINE  b 
Cl 
C2 
C3 

CE 

LINE  7 
Cl 
C2 
C3 

LINE  a 

Cl  Number  of  days  to  be  modelled  plus  one. 

LINE  9  (Dally  data) 

Cl  Julian  date 

C2  time  of  observation  (24  hour  notation) 

C3  air  temperature  (  C) 

C4  relative  humidity  (T) 

C5  cloud  cover  (teijths) 

C6  wind  speed  (ms  ) 

C7  precipitation  (mm  water) 


soil  depth  (cm) 

start  soil  temperature  at  known  depth  (  C) 
start  soli  surface  temperature  (  C) 


area  of  cell  (metres) 
critical  snow  depth  (cm) 
number  of  seconds  in  a  day 
bcoefflcient  (decimal) 
saturation  (decimal) 


slope  angle  (degrees) 
slope  aspect  (degrees) 
start  date  (Julian  calendar) 
latitude  (decimal) 


Table  14  :  Cloud  Genera  and  Cloud  Type  Indices 


Cloud  Genera 

Abbreviation 

Index 

Value 

Comments 

Cirrus 

Ci 

I 

High  clouds  composed  of  white  delicate 
filaments,  patches  of  narrow  bands, 
elements  often  curved  or  slanted  and 
smaller  than  Cs ,  never  overcast  or 
precipitating 

Clrrostratus 

Cs 

2 

High  clouds  appearing  as  whitish  veil 
usually  fibrous,  often  produces  halo 
phenomena,  thinner  than  As,  does  not 
appear  to  move,  nonprecipitating 

Altocumulus 

Ac 

3 

Hldlevel  clouds,  patches,  usually 
broken,  lee  wave  clouds,  elements 
smaller  chan  Sc,  nonprecipitating 

A1  tostratus 

As 

Midlevel  grey  sheet  or  layer  of 
striated,  fibrous  or  uniform 
appearance,  large  horizontal  extent; 
thicker  chan  Cs ,  thinnei  than  Ns, 
precipitation  general -.y  light  and 
continuous  (if  any) 

Stratocuraulus 

Sc 

5 

Grey  and/or  whitish  layer  or  patch, 
nearly  always  has  dark  spots  and  is 
nonfibrous;  elements  larger  chan  AC, 
nonprecipl eating 

Stratus 

St 

6 

Grey  rather  uniform  base,  patches 
ragged  if  present,  precipitation 
unusual  but  light  and  contnuous  if 
present,  lower  and  more  uniform  chan 

Sc,  less  dense  and  less  "wet"  than  Ns 

Nlmbos  tratus 

Ns 

7 

Grey  often  dark,  diffuse,  large 
horizontal  and  vertical  extent,  thicker 
than  As,  more  uniform  than  Sc,  often 
precipitating,  precipitation  continuous 

Fog 

FC 

8 

*  Cloud  genera;  Cumulus  (Cu),  Cirrocumulus  (Cc),  and  Cumulonimbus  (Cb) 
are  not  treated  here.  At  low  cloud  covers  (  0.3),  Cu  and  Cc  may  be 
approximated  with  Ac. 
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C4  Latitude:  The  latitude  of  the  cell  expressed  as  a  decimal, 

LINE  6 

Cl  Cell  area:  This  is  calculated  by  a  manual  subdivision  of  the 

catchment.  Subdivision  criteria  are  slope  aspect,  vegetation  cover 
and  altitude. 

C2  Critical  snowdepth:  This  is  the  snowdepth,  5  cm,  which  demarks  the 
presence  and  absence  of  snow,  as  discussed  in  section  3.2. 

C3  DELTS:  The  number  of  seconds  in  a  day,  used  in  SNOMO  calculations. 

C4  Beta  coefficient:  Parc  of  Che  energy  budget  melt  equation 

(equation  6)  . 

C5  Saturation:  Percentage  saturation  of  the  snow  or  soil  expressed  as 
a  decimal . 

LINE  7 

Cl  Soil  depth  :  This  is  either  measured  In  the  field  or  taken  from 
soil  maps. 

C2  Soil  and  surface  soil  temperature:  The  soil  temperature  at 

&  either  the  base  of  the  profile  or  at  a  known  depth,  which  is 

C3  then  designated  as  the  base  of  the  soil  profile.  Soil  temperatures 

are  measured  in  the  field  (soil  probe  or  frost  tube).  Soil 
temperatures  at  depth  are  required  only  if  SNOMO  is  initiated  from 
a  bare  ground  situation.  Soil  temperature  could  be  roughly 
estimated  from  the  air  temperature  and  preceding  air  temperatures 
if  necessary.  SNOMO  operates  on  the  extrapolation  of  the  daily 
Input  data  on  a  dally  basis.  Thus,  for  Che  extrapolation  to  occur 
the  current  day's  and  the  next  day's  data  has  to  be  input. 
Therefore,  for  the  number  of  days  of  dally  data  input  (n)  there  are 
always  n-1  days  modelled  by  SNOMO. 
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LINE  8 

Cl  n,  the  number  of  days  to  be  modelled  plus  one.  The  additional  day 
allows  for  the  data  extrapolation  routines  in  SNOMO. 

LINE  9  -  daily  data 

Cl  Julian  date 

C2  Time  of  observation.  This  is  input  in  the  24-hour  clock  format. 

Time  of  observation  is  used  as  this  allows  for  flexibility  in  time 
input,  useful  for  remote  manual  sites. 

C3  Air  temperature. 

C4  Relative  humidity. 

C5  Cloud  cover  is  measured  in  tenths  and  input  as  a  decimal.  If 

observation  data  is  unavailable  then,  as  with  cloud  type,  it  is 
estimated  from  synoptic  conditions. 

C6  Precipitation  is  that  which  is  collected  In  a  standard  snow/rain 

gauge,  l.e.  it  Is  measured  in  millimetres  of  water  and  not  snow 
depth.  SNOMO  distinguishes  between  rain  and  snow  precipitation 
after  the  precipitation  data  Is  input. 

Data  is  also  Included  in  the  SNOMO  code  in  DATA  statements.  This  is 

physical  constants  and  computational  control  data: 

(a)  Density  of  water,  1000  kgm 

(b)  Heat  capacity  of  water,  4.21  x  10^  Jm  ^  k  *. 

6  ”3  ■“  1 

(c)  Heat  capacity  of  snow,  0.21  x  10  Jm  k 

(d)  Latent  heat  of  fusion,  0.334  x  10^  Jkg 

(e)  TOTTIM:  Total  number  of  24  hr  repetitions  used  in  solving  the  TSTM 
heat  flow  calculation.  This  is  set  to  1. 

(f)  TFRQ:  Time  step  In  minutes,  used  in  solving  the  TSTM  heat  flow 
calculation,  set  at  5  minutes. 
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(8) 


TPRNT:  Output  time  frequency  In  minutes.  This  refers  to  the 
TSTM  outputs  of  surface  temperature  (°C),  greybody  radiation  (Wm”^) 

solar  insolation  (WM  ),  surface  absorption  (Wta  ),  atmospheric 
—2 

infra-red  emission  (Wm  ),  sensible  heat  (Mm  ")  and  latent  heat 
(Wm"2). 

Appendix  I  summarises  the  logic  and  operational  structure  followed  by 
SNOMO.  The  logic  and  operational  structure  followed  by  TSTM  is  contained 
in  Appendix  II. 

3. 4  Calculation  of  the  snow/soil  energy  budget 

SNOMO  is  being  developed  to  model  the  snow  system  and,  therefore,  unless 
specified,  the  remainder  of  this  section  is  concerned  with  the  snow 
system  rather  than  the  soil  system.  To  avoid  confusion  and  repetition 
TSTM(l)  will  be  used  to  indicate  the  original  TSTM  program  (Balick  et  al . 
1981a)  and  TSTM(2)  to  indicate  the  converted  TSTM  subroutine  used  in 
SNOMO.  When  TSTM  Is  used.  It  will  refer  to  TSfM(l)  and  TSTM(2). 

The  snow  energy  budget  is  calculated  using  the  subroutine  TSTM.  T>><s,  as 
stated  previously,  was  converted  from  the  Terrain  Surface  temperature 
Model  (Balick  et  al.  1981a)  with  the  help  of  Mr  Randy  Scoggins  (WES). 

Figure  37  shows  the  concept  for  TSTM.  TSTM  predicts  surface  temperatures 
for  a  multilayered  (1-6  layers)  system  by  determining  energy  transfer  in 
and  out  of  the  system.  The  model  assumes  chat  the  layers  and  Che 
environment  above  them  are  horizontally  uniform,  i.e.  the  most 
significant  heat  fluxes  ar  e  vertical.  Therefore,  the  temperature  T 
estimates  result  from  solving  the  one-dimensional  heat  equation. 

£)  =  yr(i, 

Subject  to  the  boundary  conditions. 

> 

^  =  0  oJb  €  »  0 

i.=) 


— •  (16) 


(17) 


and 


(18) 


L=l 

Where  the  observable  surface  Is  z  =  b,  lower  surface  is  z  =  B,  (z)  is 
Che  diffusivlty  and  both  and  8^^.,  1  =  1,2,..,  n  denote  heat  fluxes 
and  at  time  t.  Reliability  of  the  results  depends  upon  the  extent  to 
which  Che  thermal  characteristics  in  each  of  the  layers  can  be 
approximated  by  constant  values  and  Is  strongly  dependent  upon  the 
approximation  of  b-i  .  1  =  l,2,...,n,  taking  place  at  the  surface  which  is 
exposed  to  environmental  heat  fluxes.  TSTM  calculates  and  outputs  at  any 
desired  time  interval  within  24  hours  (in  TSTM(2)  this  has  been  changed 
to  any  time  interval  between  one  observation  time  and  the  next)  Che 

surface  temperature  (°C),  greybody  radiation  Wm  ^),  solar 

—2 

insolation  (K^,  Wm  ')  ,  surface  absorption  (KJ/  -  K'^  ,  Wm  )  ,  atmospheric 

—2  “ 

Infra-red  emission  (L'l',  Wm  ),  sensible  heat  loss  (Aq^  ,  Wm  *■)  and  latent 
heat  loss  (AQ^i  Wm  )  of  Che  snowpack. 

TSTM(l)  was  designed  Co  operate  on  a  larger  data  base  than  TSTM(2). 
TSTM(l)  was  designed  to  allow  the  input  of  solar  insolation  but  with  the 
calculation  of  this  if  unavailable.  The  input  of  measured  solar 
insolation  values  was  preferred.  TSTM(2)  has  been  modified  so  that  only 
calculation  of  solar  Insolation  occurs  but  with  little  modification 
TSTM(2)  could  revert  back  to  the  input  of  solar  insolation  values.  This 
could  possibly  be  incorporated  into  an  initial  user  menu. 

3.4.1  Calculation  of  incoming  solar  radiation, 

TSTM  calculates  only  the  direct  incoming  solar  radiation.  The  solar 
radiation  incident  at  the  cop  of  the  atmosphere,  the  solar  constant,  S^, 
is  depleted  by  the  atmosphere  depending  upon  the  length  of  its  path  and 
the  vertical  transmissivity  of  the  air.  The  path  length  (or  optical  air 
mass  number,  M)  is  the  ratio  of  the  slant  path,  taken  by  the  beam  to  the 
zenith  distance,  so  that  M  -  secz  •  l/cosz.  The  effects  of  altitude  are 
allowed  for.  Atmospheric  transmissivity  depends  upon  the  concentrations 
of  gases,  droplets  and  particles  in  the  atmosphere. 
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Therefore,  S,  the  direct  beam  solar  radiation  on  a  horizontal  surface 
(excluding  the  effect  of  clouds  and  optical  air  mass  number)  is: 


S=  a)](0-2il)S.  CoS  t 

C0S2r 


where 

surface  albedo 

Mugge-Moller  absorption  function,  equal  to  0.271  (u*  secz)^’^*^^ 
effective  water  vapour  content  of  atmosphere 
(p*54^^5^mount  of  solar  radiation  of  wavelength  greater  than  0.9  urn 
solar  radiation  incident  on  top  of  the  atmosphere 
2r  zenith  angle  of  the  sun  as  a  function  of  time  of  the  day  and  time 
of  the  year 

0^^  atmospheric  albedo  for  Rayleigh  scattering,  equal  to  0,085-0.247 

^°Sl0 

surface  pressure 
lOOOmb 

5^  area  average  ground  albedo 

solar  radiation  of  wavelength  less  than  0.9  um 

The  effective  water  vapour  content,  f* ,  is  the  total  precipitable  water 
excluding  clouds.  Precipitable  water  is  estimated  from  surface  air 
temperature  and  relative  humidity 

=  exp  1^0 -0^0^ 4- 7^  (20) 

where 

Tdl  dew  point  temperature  (°C) 

-  0.02290  April  -  June 
0.02023  July  -  March 
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The  effect  upon  S  of  cloudy  skies  is  calculated  in  two  stages.  First,  a 
cloud  adjustment,  factor,  CA,  is  calculated  which  is  dependent  upon  cloud 
type.  This  is  then  related  to  cloud  cover  to  result  in  S^. 

CA  =  (a/94.4)  X  exp  [-m  X  (  b-0.059)j  (21) 

where 

a  and  b  =  empirical  coefficients  dependent  upon  cloud  type  (table  15) 
m  =  optical  air  mass  number 


S 

c 


-{Ls- 


(S  X  CA)}  X 


cc 


where 

CC  visual  cloud  cover  in  tenths 

is  then  modified  for  slope  and  aspect  by  the  calculation  of 
a  slope  factor,  SF 

k1  =  S  X  SF  (23) 

^  c 


spherical  trigonometry  gives  the  following  relationships,  figure  33. 


COSir  Sinjjf  +  COSJ^COsS  OOsU  =  siki 

(si/i^COSjZ^  -  COs5sM|Z^COsk)/SlVl  2 


=  3()0°^  (si/iS’cos^  -  CosS'stVv^  (r»sk.)/sM  -2. ,  012. 
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Table  15  :  Coefficients  used  in  SNOMO  energy  budget  calculations 


Cloud  type 

Coefficient 

a 

b 

CIR 

Cirrus 

82.2 

0.079 

0.04 

Cirrostratus 

87.1 

0.148 

0.08 

A1  cocumulus 

52.5 

0.112 

0.  17 

A1  tostratus 

39.0 

0.063 

0.20 

Stratocumulus 

34.7 

j.  104 

0.22 

Stratus 

23.8 

0.159 

0.24 

Nimbostratus 

11.2 

-0.  167 

0.24 

Fog 

15.4 

'•'.028 

0.25 
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where 

i 

£ 

S 

A 

k 


apparent  solar  time 
solar  zenith  angle 
latitude  of  location 
solar  declination 
solar  azimuth  angle 
hour  angle 


For  a  slope,  the  solar  radiation  Incident  on  the  slope,  S,  can  be 
calculated  using  spherical  trigonometry,  fl,.  ire  39. 


COS@=  C0S2r4  StVl^  SiVl?r  cos(JZ.-il) 


where 


4 

COS0 


slope  angle 
slope  azimuth  angle 
slope  factor 


3.4.2  Calculation  of  reflected  solar  radiation, 

TSTM  calculates  the  amount  of  direct  Incoming  solar  radiation, 
is  absorbed  at  the  snow  surface,  K  : 

3  D 


where 


oC 

0-a) 


surface  albedo 
absorptivity  of 


the 


surface 


(27) 


that 


(28) 
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The  reflected  solar  radiation,  can  therefore  be  calculated; 


U'ts  (29) 

or  alternatively, 

Kf  =  K'l'oC 

(30) 

3.4.3  Calculation  of  incoming  longwave  radiation  (ataospherle  Infra-red 
emission)  Li 

The  Brunt  equation  is  used  to  calculate  atmospheric  IR  radiation  on  the 

surface ,  li  : 

o 

II- £(rT<,‘^[t+t(eDJ  <„, 

where 

^  emlsslvlty,  assumed  equal  to  1 

(T  Stephan-Bol tzmann  constant  (0.813  x  10  Cal  cm  ”  min  *  °rC  '*) 

*f^  shelter  air  temperature,  °K. 

S.^  water  vapour  pressure,  mb 
b  y  C.  empirical  constants  =  C  =  0.61,  h  =  0.05 

The  value  of  is  obtained  from  Teten's  equation: 

e^=  ^.|02  x0xp(.AKT^)/(rct  ■H2^3-IS'-g) 


where 


A 

3 


relative  humidity  (decimal) 

17.269 

35.86 


?C 


Clouds  also  contribute  to  L'^  : 


Where 

CI^ 


coefficient  dependent  upon  cloud  type  (see  table  15). 


(33) 


3.4.4  Calculation  of  reflected  longwave  (greybody  radiation,  if 


where 


emissivity  of  ground 

current  surface  temperature  as  predicted  by  the  model 


(34) 


3.4.5  Calculation  of  sensible  heat  flux,  ^  Qj, 

^  ^  SCP 


where 

r  l.I75(l-15Rl)°*^^ 

Ri  ^  0 

SCFJ  (1.5R1)^ 

0<R1  ^0.2 

i  0 

Ri  >  0.2 

(35) 


^v. 


M 

bt 


Q 


■(?c 


air  density 

specific  heat  of  dry  air  at  constant  pressure 
von  Karman's  constant  (0.40) 
observation  height 

partial  derivative  of  potential  temperature 
with  respect  to  height  z 

partial  derivative  of  wind  speed  with  respect  to  height 
Richardson  number 
potential  temperature 


z 


4^ 
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Potential  temperature , W  ,  is  defined  by: 


where 

air  temperature 
P  air  pressure 

The  Richardson  number,  R1 ,  is  defined  by: 


where 

^  acceleration  due  to  gravity 

^  average  potential  temperature  between  the  surface 
and  height  z 


(36) 


(37) 


In  TSTM  and  ^'/^^iare  approximated  by  first  order  differences  and  It 

is  assumed  that  the  air  temperature  at  the  surface  equals  temperature  of 
the  surface  and  that  wind  velocity  at  the  surface  is  zero. 


3.4.6  Calculation  of  latent  heat  flux.AQ,, 

^  SCF  (38) 

where 

L.  latent  heat  of  evaporation,  597.3  cal  g  ^ 

^  specific  humidity 

w  decimal  relative  saturation  of  the  top  surface 
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3.4.7  Calculation  of  surface  temperature  and  solution  of 


heat  flow  equations 

At  present,  TSTM(2)  assumes  a  constant  temperature  at  the  bottom 
boundary  (the  snow/soll  Interface).  TSTM(l)  allows  one  of  the  following 
three  options: 

(a)  Option  1  A  constant  temperature  at  the  bottom  boundary. 

(b)  Option  2:  A  constant  heat  flux  at  the  bottom  boundary. 

(c)  Option  3:  A  constant  heat  flux  at  the  bottom  boundary  and  an 

additional  constant  temperature  radiating  surface  below  the  bottom 
boundary. 

Option  3  requires  additional  Input  data:  I'l)  bottom  boundary  thermal  IR 
emlsslvlty;  (11)  bottom  boundary  geometric  shape  factor;  (111)  under 
surface  thermal  IR  emlsslvlty;  (tv)  under  surface  geometric  shape 
factor;  and  (v)  under  surface  temperature.  The  bottom  boundary 
condition  Is  kept  consCanc  In  time  regardless  of  the  option  chosen. 

The  complicated  nonlinear  boundary  conditions  require  that  the  heat 
conduction  equation 

~ir~  — 

be  solved  numerically.  In  this  equation,  oC  Is  the  diffusivlty.  Each 
layer  Is  assumed  to  be  homogeneous,  and  it  is  assumed  that  the  thermal 
characteristics  can  be  taken  to  be  constant;  specifically,  the  thermal 
conductivity  and  diffusivlty  for  each  layer  are  assumed  to  be  a  constant. 

3.4.8  Solution  within  a  layer 

Within  each  layer,  an  explicit  scheme  is  employed  to  solve  the 
one-dlmenslonal  heat  equation.  In  particular,  given  the  temperature 
profile  at  time  t  ,  the  temperature  at  time  t  +At  at  the  node  z  Is 
given  by 


(39) 


where  A  t  Is  the  time  increment  andAz  denotes  the  spatial  increment. 

It  should  be  noted  that  numerical  stability  requires 

The  problem  of  numerical  stability  is  critical  for  thin  highly  conductive 
layers . 

3.4.9  Solution  at  the  interface  of  two  layers 

The  following  derivation  of  an  explicit  finite  difference  scheme  to 
handle  the  interface  between  layers  is  a  modification  of  that  presented 
in  Carnahan,  Luther,  and  Wilkes  (1964).  The  derivation  assumes  perfect 
thermal  contact  at  the  interface,  I.",  continuity  of  the  IieaL  flux  and 
temperatures  at  the  Interfaces. 

Let  layer  1  have  thermal  conductivity,  k^,  and  dif  fuslvity  ,  and  layer 

2  have  thermal  conductivity  and  diffusivity,  k^,  and  <lL2>  respectively. 


Knowing  the  temperature  Tj^_p  T^^,  and  T^^^j  at  the  node  points  1  -  1,  i, 
and  1+1,  the  problem  is  to  calculate  the  new  temperature  Tj^  at  the 
Interface.  Employing  the  truncated  Taylor  series,  Tj^_j  is  approximated 
by 


(41) 


94 


where  il  denotes  the  partial  derivative  In  layer  1  at  the  Interface. 
Thus,  ^ 


Also,  the  first  order  approximation  to  Is  given  by 


Since  one  obtains 


'  v^ki 


(42) 


(43) 


(44) 


In  a  similar  fashion  for  layer  2,  one  obtains 


..  [-rl 

rj  f 

J  4., 
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Continuity  of  the  heat  flux  implies  that  Equation  46  equals  Equation  47; 
thus,  after  simplification,  the  final  equation  used  to  calculate  is: 


,  1: 

I  k  ~1r  lii  f  _/^ 


(48) 


3.4.10  Upper  boundary 


The  new  or  updated  value  of  the  surface  temperature  T(t  +  t,0)  is 

calculated  by  solving  the  surface  heat  balance  equation  for  the  surface 
temperature,  T^.  The  beat  balance  equation  is 


K*  +*  +■  +■  ^ 


(49) 


whe  re  denotes  the  heat  flux  Into  the  surface;  i,e,  * /^2r) 

and  is  approximated  by  |j  /Al)  where  k.  denotes  the  conductivity  of 

the  surface  layer  and  T^  denotes  the  temperature  at  the  present  time  for 
the  first  node  point  below  the  surface.  Letting 
the  heat  balance  equation  becomes 


-£6'l^f-hk 


(50) 


or  upon  rewriting 


£j6'Ai  mi  its 
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The  function  F  is  defined  by 


F(rj)=T^’'+ 


iWiro, 


It  is  seen  that  the  updated  surface  temperature  is  a  root  of  F.  The 
Newton  Faphson  algorithm  has  been  employed  to  locate  a  value  of  T^  such 
that  F  vanishes.  In  employing  the  Sewton-Raphson  scheme,  the  derivative 
of  F  with  respect  to  T^  is  needed: 


r  4-Tj  -i- 

iTo,  J  ^6'Ai 


Numerical  considerations  have  resulted  in  the  approximation  of  dD/dT^ 
by  the  following  expression 


(ib  _  (  t>h/  -  ^0^ 


where  D„  is  the  value  of  D  using  the  latest  estimate  of  T  ,  D,  is  the 
N  °  g  0 

value  of  D  obtained  by  using  the  previous  estimate  of  T^ ,  and  T  denotes 
the  change  in  temperature.  The  starting  value  for  the  Newton-Raphson 
scheme  is  taken  to  be  the  surface  temperature  at  the  previous  time  step. 
It  appears  that  three  to  five  iterations  yield  satisfactory  convergence 
to  the  new  surface  temperature. 


3.4.11  Bottom  boundary 

The  bottom  boundary  condition  is  the  heat  flux  through  the  bottom  of  the 
lowest  layer  and  can  be  specified  with  one  of  three  options.  The 
requirement  of  a  constant  temperature  results  in  a  straightforward 
boundary  condition. 


For  options  2  and  3,  it  is  required  that  the  following  equation  be 
satisfied: 
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where  denotes  the  radiative  energy  loss  through  the  bottom  boundary,  G 
denotes  the  heat  flux  into  the  lower  surface  and  is  given  by 
where  k  denotes  the  conductivity  of  the  bottom  layer,  R^^  denotes  the 
radiative  energy  from  the  constant  temperature  radiating  surface  below 
the  bottom  boundary,  and  D  is  the  constant  heat  flux  at  the  bottom 
boundary.  G  is  approximated  by  where  Tg  is  the 

temperature  of  the  bottom  surface  and  Tj  is  the  temperature  at  the  first 
node  point  above  the  bottom  surface. 


The  following  equation  results  from  substituting  the  appropriate  energy 
components  into  Equation 

where  Cg  denotes  the  bottom  boundary  thermal  IR  emlssivity,  b^^g  denotes 
the  bottom  geometric  shape  factor,  is  the  under  surface  thermal  IR 
emlssivity,  bj^j^  is  the  under  surface  geometric  shape  factor,  and  Tj^ 
denotes  the  under  surface  temperature.  Equation  ^4  is  solved  by 
employing  the  Newton-Raphson  iterative  scheme. 


3.5  Raln-on-snow 


The  heat  transferred  to  the  snow  by  rainwater  is  the  difference  between 
its  energy  content  before  falling  on  the  snow  and  its  energy  content  on 
reaching  thermal  equilibrium  within  the  pack.  Depending  upon  the  thermal 
state  of  the  pack,  the  energy  produced  by  the  rain  will  heat  up  the  pack 
or  melt  the  pack. 

1)  Rainfall  on  a  melting  pack  where  the  rainwater  does  not  freeze 

Op-/2Cp(rr-roP/ia)0  ,37, 
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where 


—2  - 1 

energy  supplied  to  the  pack  by  rain,  kjm  day 
density  of  water 

heat  capacity  of  water,  kJkg”^*^C  ^ 
temperature  of  the  rain,  °C 
temperature  of  the  snow,  °C 
precipitation  rate,  mm  day  ' 


(from  Male  and  Gray  1981) 


In  SNOMO  T  is  taken  to  be  the  same  as  the  air  temperature  T..  In  a 
r  A 

melting  pack,  the  pack  should  be  isothermal  at  0  C,  therefore  the  working 
equation  for  SNOMO  is: 


(5p=  P)/lOOO 


li) 


where 

P 

4 

V 

4 

Tw 

Pc 

2> 


Rainfall  on  a  pack  with  a  temperature  below  0  C,  where  the  rainfall 
freezes  and  releases  its  latent  heat  of  fusion  resulting  in  a 
warming  of  the  pack.  The  rate  of  change  in  snow  temperature,  irjit 
Is  calculated  from  Male  and  Granger  (1978): 


-Cp( 

ctb  Cp,p.i> 


precipitation  rate,  ram  day 

_3 

density  of  water,  kgm 

heat  capacity  of  water,  kjkg  *  °C  * 

temperature  of  rain,  taken  to  equal  air  temperature,  °C 

latent  heat  of  fusion,  kJkg  * 

heat  capacity  of  snow,  kjkg  ^  °C  * 

average  snow  temperature ,  °C 

snow  density,  kgm 

depth,  m. 
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3.6  Calculation  of  meltrate 

The  calculation  of  meltrate  occurs  within  the  two  subroutines  OIELT  and 
CIMELT.  CMELT  and  CIMELT  are  Identical  except  that  CMELT  uses  dally 
energy  flux  totals  and  CIMELT  uses  average  dally  energy  fluxes  In  the 
calculation  of  the  energy  available  for  mell.^Q^j*  CMELT  and  CIMELT  are 
used  for  comparative  purposes  In  the  development  of  SNOMO.  The  meltrate 
calculation  Involves  five  stages: 

1)  Calculation  of  snow  density,  P 

l  sn 

If  the  snowpack  is  single  layered  then  the  snow  density  is  the 
value  in  THKMXl  (1,6).  If  the  snowpack  is  two-layered  then  the 
average  snow  density  is  calculated  allowing  for  differing  layer 
thicknesses . 

11)  Calculation  of  the  average  temperature  of  the  snowpack,  T„ 

M 

This  Is  calculated  using  the  daily  average  surface  and  base 
temperatures  of  the  snowpack. 

Ill)  Calculation  of  the  snowpack  Internal  energy  change,  dU/dt 


di#  --  )>[p.  Cf>. 


aju/air  snowpack  internal  energy  change,  kJra  “  liav 
^  snowdepth,  m 
^  density,  kgm  ^ 

Cp  specific  heat,  kjkg  ^  °C  ' 

Mean  snow  temperature,  °C 


and  1,  1,  V  refer  to  Ice,  liquid  and  vapour 


(from  Male  and  Gray,  1980 
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If  It  Is  assumed  that  the  humidity  of  the  air  in  the  snowpack  Is 
100%  then  the  contribution  to  dU/dt  by  the  vapour  phase  term  Is 
negligible  and  is  Ignored  in  the  calculation  of  du/dt.  During 
non-melt  periods,  there  is  no  liquid  component  to  Che  pack  and, 
therefore,  terms  and  are  not  required.  The  working  equation 
for  a  melting  pack  is  therefore; 

dU/db  --  ( b n 

Iv)  Calculation  of  the  energy  available  for  melt.Ag 

M 

This  calculation  is  based  on  equation  (5)  and  o,.curs  in  i  stages. 
Dally  totals  or  averages  are  used. 

a)  Q*  »  K-V-  Kf+  L'l/-  Lt  (62) 

where 

0*  net  radiation 
K't  incoming  shortwave  radiation 
Kt  reflected  shortwave  radiation 
L'V  incoming  longwave  radiation 
ft  reflected  longwave  radiation 

_  t 

Units;  Wm 

K^r,  and  L'Tare  calculated  by  TSTM  and  KVis  solved  from  the 
TSTM  term  XABSOR 


b)  Using  equation  (5): 


ENAVM  -  Q*  +  Q,_  +  Q  +  Q 
h  e  g 

where 

ENAVM  energy  available  for  meit,  intermediate  term 
sensible  heat  flux,  ATERM 
Qg  latent  heat  flux,  BTERM 

Qg  ground  heat  flux,  set  to  zero  at  present 

Units;  Wra“^ 


(63) 
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— —2  “I 

ENAVM  is  converted  from  Wm  “  to  Jm  day  by  multiplying  by 
8.64  X  10^. 


c)  For  the  completion  of  equation  (5): 


AQv, 


ENAVM  +  Qp  ■  dU/dt 


(64) 


where 

A  -2-1 

QQy,  energy  available  for  melt,  Jm  day 

^  -2-1 
Q  energy  Introduced  to  the  pack  by  rain,  Jm  day 


v)  Calculation  of  the  meltrate,  M 


The  meltrate  is  calculated  from  equation  (6)  in  both  millimetres  of 
water  equivalent  and  centimetres  of  snow. 


h 


AQ„/(/f^jl.looo 


■1 


where 

Njjj  meltrate,  mm  water  day 

latent  heat  of  fusion,  Jkg  * 
density  of  water,  kgm  ^ 


(65) 
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where 

meltrate,  cm  snow  day-1 
-3 


Sh 


^  density  of  snow,  kgm 


(66) 


3.7  Snowfall  and  compaction  changes  to  pack  characteristics 

This  section  of  SNOMO  models  the  effects  of  snowfall,  or  absence  of 
snowfall,  on  the  snowpack.  SNOMO  models  the  effects  of  three  snowfall 
Intervals  : 
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i)  No  snowfall  for  3  days  or  longer. 

If  3  days  after  snowfall,  no  snow  has  fallen  on  the  fourth  day  any 
'new'  snow  in  the  pack  becomes  'old'  snow. 

11)  Snowfall  for  2  consecutive  days. 

Snow  falling  on  any  2  consecutive  days  is  considered  to  belong  to 
the  same  snow  event.  Therefore,  snowfall  on  day  2  is  simply  added 
to  the  'new'  snow  total  that  fell  on  day  1. 

ill)  Time  Interval  of  1  or  2  days  between  snowfalls. 

This  is  demonstrated  by  the  situation  when  there  is  a  snowfall  o.i 
'day  1,  no  snowfall  on  day  2,  and  snowfall  on  day  3.  Here  the  'new' 
snow  that  fell  on  day  1  is  converted  to  'old'snow  on  day  3.  The 
two  snowfalls  are  considered  as  separate  events. 

Each  change  in  the  snowpack  layering  results  in  a  change  in  the 
physical  characteristics  of  the  snowpack,  due  to  the  substitution 
of  VNEWSN  and  VOLDSN  values  in  to  THKMXl  or  THKMX2.  Therefore,  the 
pack  characteristics  of  thermal  diffuslvity,  heat  conductivity, 
emlsslvlty,  albedo,  and  density  are  modelled  as  changing  with  time 
and  snowfall  (frequency  and  amount).  When  'new'  snow  is  converted 
to  'old'  snow  the  resultant  'old'  snow  layer  is  thinner  than  the 
original  'new'  snow  layer.  This  is  because  of  the  higher  density 

-3 

of  the  'old'  snow  (0.48  gm  )  compared  to  the  'new'  snow  (0.10 
gtn”^) 

SNOMO  at  present  uses  the  extreme  values  for  the  physical 
characteristics  of  the  pack,  l.e.  the  values  for  freshly  fallen 
snow  and  very  old,  dirty  snow.  When  the  'new'  to  'old'  snow  or 
'old'  to  'new'  snow  conversion  occurs,  this  crer.tes  large  changes 
in  the  snowpack  characteristics  overnight.  A  situation  not 
replicated,  except  for  albedo,  in  reality.  However,  this  method  is 
sufficient  for  SNOMO  at  present.  It  would  not  be  very  difficult 
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to  Introduce  a  seasonal  variation  In  certain  values,  e.g.  density, 
and  to  lessen  the  effect  of  overnight  changes, 

3.8  Melt  changes  to  pack  characteristics 

This  section  of  SNOMO  models  the  effect  of  melt  on  the  snowpack. 

Meltrate  in  cm  snow  day”''  and  mm  water  day”'  (where  day  refers  to  a  SNOMO 
operational  day)  is  calculated  by  SNOMO.  The  depth  of  snow  subtracted 
from  the  pack  if  melt  occurs  is  measured  in  centimetres  (i,e.  calculated 
using  snow  density,  equation  66)  and  operates  over  one  operational  day. 

A  threshold  melt  depth  of  1  cm  is  used  below  which  no  melt  is  said  to 
have  occurred.  The  subtraction  of  any  melt  is  dependent  upon  the  layer 
structure  of  the  pack  with  substractlon  always  from  the  surface  of  the 
pack  downwards.  Therefore  if  melt  occurs  in: 

i)  a  1-layer  pack 

-  SNOMO  stops  if  the  melt  depth  is  greater  than  the  pack  depth 
subtraction  of  the  melt  depth  occurs  and  the  pack  remains 
single-layered 

li)  a  2-layered  pack 

SNOMO  stops  if  the  melt  depth  is  greater  than  the  pack  depth 

-  if  the  melt  depth  is  greater  than  the  thickness  of  the  top 
layer,  then  the  melt  depth  is  subtracted  from  the  surface 
resulting  in  a  single-layered  pack 

if  the  melt  depth  is  less  than  the  thickness  of  the  top  layer 
then  subtraction  occurs  and  the  pack  remains  double-layered 

3.9  Output 

At  present,  the  output  from  SNOMO  reflects  the  current  phase  of 
development.  Once  fully  operational,  the  output  of  SNOMO  will  be 
manipulated  so  that  it  is  menu  driven,  or  at  least  offers  the  operator  a 
choice  of  outputs.  The  output  could  be  produced  in  graphical  form. 
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IV  CURRENT  RESEARCH 
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IV.  CURRENT  RESEARCH 


4. 1  Areas  of  current  research 

There  are  four  areas  of  current  research: 

layering  and  the  manipulation  of  snowpack  characteristics  have  been 
Introduced  to  SNOMO  and  are  now  operational.  Initial  comparisons 
between  calculated  radiation  results  and  measured  results  at 
Townline  (W3,  Danville,  Vt.)  are  favourable. 

SNOMO  currently  calculates  melt  depth  at  the  cell  scale  but  the 
timing  of  the  melt  requires  further  modification. 

WES  visit  (27  February  -  3rd  March  1989). 

1)  VEGIE  (Ballck  et  al.  1981b)  was  incorporated  into  SNOMO  and 
made  operational  with  the  help  of  Mr  Randy  Scoggins.  The 
inclusion  of  VEGIE  enables  the  modelling  of  the  effects  of 
vegetation  (up  to  one  metre  in  height)  on  the  snowpack. 

11)  The  interpolation  routine  used  by  TSTM2  to  calculate  air 
temperature  from  the  two  input  points  was  Improved.  The 
interpolation  now  reflects  the  diurnal  temperature  cycle. 

iil)  The  TSTM(2)  code  was  optimized  removing  any  unused  remnants 
of  TSTM(l). 

CRREL  visit  (6th  March  -  7th  April  1989). 

i)  Discussed  SNOMO  with  Dr  Rachel  Jordan  (CRREL)  who  suggested 
various  improvements. 

11)  Field  program  conducted  at  W3,  Danville,  Vermont.  The  object 
of  the  field  program  was  to  maximize  the  collection  of 
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spatial  snow  depth  and  water  equivalent  data  over  W3  during 
the  snowmelt  period.  A  series  of  trails  were  laid  across  W3 
and  snow  courses  set  at  points  on  the  trails  to  reflect  the 
vegetation,  aspect  and  altitude  differences  over  W3.  In 
addition  pre-existing  sites  and  sites  used  in  1988  were  also 
Included.  Every  Tuesday  for  4  weeks  data  from  40  snow 
courses  (usually  5  sampling  points  per  snow  course)  was 
collected.  In  addition  to  Dr  Timothy  Pangburn,  Mr  Bill 
Roberts  and  Mr  Hugh  Greenan,  extra  researchers  came  from 
CRREL.  The  data  collected  will  be  used  to  validate  SNOMO  but 
will  also  be  used  at  CRREL  (Dr  Timothy  Pangburn  and  Dr  Ike 
Mackim)  and  at  UVM  (Dr  Alan  Cassell). 
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OF  RESEARCH  PROJECT  ACHIEVEMEHTS 
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SUMMARY  OF  RESEARCH  PROJECT  ACHIEVEMENTS 


Development  of  a  phototype  distributed  snowmelt  model,  SNOMO,  based 
on  the  energy-budget  approach  to  snowmelt  calculation.  The 
development  has  Involved: 

1)  development  of  the  SNOMO  operational  logic  allowing  for  the 
manipulation  of  the  TSTM  subroutine  and  the  effect  of  melt 
and  snowfall  on  the  temporal  and  spatial  variations  of  the 
snowpack  characteristics 

li)  implementation  of  this  logic  into  a  'control'  program  and  the 
conversion  of  TSTM  into  a  subroutine  within  SNOMO 

iii)  incorporation  of  VEGIE,  an  addition  to  the  original  TSTM,  to 
SNOMO.  This  enables  the  modelling  of  the  energj-budget  of  a 
snow-covered  vegetated  surface  (the  vegetation  is  less  than 
1  m  high) 

Fieldwork  conducted  during  the  1988  and  1989  melt  seasons  has 
resulted  in  the  collection  of  spatial  snow  depth  and  water 
equivalent  data  for  use  in  the  validation  of  SNOMO.  The  data  has 
also  been  used  for  other  research  projects  at  CRREL. 
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Figure  II  ;  Simplified  flowchart  for  TSTM  subroutine, 
Balick  et  al.,  1981a.  (Sheet  1  of  5). 
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Figure  II  ;  (Sheet  2  of  5) 
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Table  I  :  Variable  definitions  :  SNOMO 


Variable 

Definition 

AIRTP 

Air  temperature,  °C 

AREA 

2 

Area  of  cell,  m 

AVABSOR 

Average  dall^  shortwave  radiation  absorbed  at  the 
surface,  Wm 

AVATERM 

-2 

Average  dally  incoming  longwave  radiation,  Wm 

AVBOTTP 

Average  dally  temperature  at  the  base  of  the  snowpack 

AVDTERM 

-2 

Average  dally  evaporative  heat  flux,  Wm 

AVGBR 

-2 

Average  dally  reflected  longwave  radiation,  Wm 

AVHTERM 

_2 

Average  daily  sensible  heat  flux,  Wm 

AVSOL 

-2 

Average  dally  incoming  shortwave  radiation,  Wm 

AVSTP 

Average  snow  temperature,  °C 

COEFB 

Substitute  constant  for  du/dt  can  be  used  in  melt 
calculation 

CRITDP 

Critical  snowdepth,  5cm 

CUMTM 

Cumulative  time,  i.e.  length  of  computational  day 

DAYl 

Julian  date  of  day  of  SNOMO  Initiation 

DDATE 

Julian  date  of  day  being  modelled 

DENSN 

-3 

Snow  density,  gm 

DENW 

Density  of  water,  kga  ^ 

DEPMXl 

Depth  matrix  for  1  layer  pack 

DEPMX2 

Depth  matrix  for  2  layered  pack 

DEPTH 

Snowpack  depth,  cm. 

DEPTH 1 

Used  in  calculation  of  DEPTH  when  converting  new  snow 
to  old  snow,  cm. 

/ cont ... 

(Sheet  1  of  3) 


Iir-2 


DUDT 

ENAVM 

ENAVMl 

GTERM 

HCAPSN 

HCAPW 

K1 

K2 

K3 

KIO 

K25 

KELVIN 

CONTROL 

LHEAF 

MRSN 


MRWE 


MRSNl 

MRWEl 

NETRAD 

NKONT 

PPTN 

PTERM 


du/dt,  part  of  Che  energy  budget  equation 

Energy  available  for  melt  Q^.  Calculated  using 
energy  totals,  Wm 

As  ENAVM  but  calculated  using  energy  averages 

Ground  heat  flux,  Wm 

Heat  capacity  of  snow,  kj  kg  ^  °C  ' 

Heat  capacity  of  water,  kJ  kg-'  °C  ' 

Count  of  the  number  of  days  modelled,  i  .e .  no.  times 
DO  5001  loop  called 

k2-l  indicates  Chat  SNOMO  modelled  soil  in  Che 
previous  day  or  is  modelling  the  first  day 

k3-l  indicates  that  soil  was  modelled  in  the  previous 
day 

Count  of  Che  number  of  internal  iterations  of  TSTM  in 
order  to  calculate  daily  averages  of  TSTM  calculated 
values. 

Count  of  Che  days  '’lapsed  since  last  snowfall 

Air  temperature  converted  to  degrees 

Counts  of  the  number  of  days  modelled  and  terminates 
SNOMO  when  necessary 

Latent  heat  of  fusion,  MJ  kg  ' 

Meltrate  in  cm  of  snow  per  'day' ,  calculated  using 
ENAVM 

Meltrate  in  ram  of  water  per  'day,  calculated  using 
ENAVMl. 

As  MRSN,  but  using  ENAVMl 

As  MREW,  but  using  ENAVMl 
-2 

Net  radiation,  Wm 

Used  with  KONTROL  to  terminate  SNOMO 
Precipitation  amount,  MM  water 

-2 

Heat  input  to  pack  by  rain-on-meltlng  snow,  Wm 

/ cont . . . 
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RH 

RSTPC 

SNPPTN 

SNVOL 

SOILDP 

SOILTP 

SSOILTP 

TABSOR 

TATERM 

TSTERM 

TGBR 

THKMXl 

THKMXa 

THTERM 

TSOL 

VNEWSN 

VOLDSN 

VSOIL 

WEI 


Relative  humidity,  X 

Rate  of  surface  temperature  change  due  to  rainfall  on 
a  frozen  pack,  °C 

Snow  precipitation  amount,  cm  snow, 

AREA  *  DEPTH,  m^ 

Soil  depth  of  base  of  soil  profile  of  depth  taken  as 
base  of  profile,  cm. 

Soil  temperature  at  depth,  taken  at  the  SOILDP,  °C 
Surface  soil  temperature,  °C 

Total  dally_|hortwave  radiation  absorbed  at  the 
surface,  Wm 

-2 

Total  daily  incoming  longwave  radiation,  Wm 

-2 

Total  daily  evaporative  heat  flux,  Wm 

-2 

Total  daily  reflected  longwave  radiation,  Wm 

Thickness  matrix  for  a  1  layer  pack 

Thickness  matrlc  for  a  2  layer  pack 

-2 

Total  daily  sensible  heat  flux,  Wm 

-2 

Total  dally  incoming  shortwave  radiation,  Wm 

Vector  holding  new  snow  values  of  thermal  diffusivity, 
heat  conductivity  emissivity,  albedo  and  density 

Vector  holding  old  snow  values  of  thermal  diffusivity, 
heat  conductivity,  emissivity,  albedo  and  density 

Vector  holding  soil  values  of  thermal  dlffusivltv, 
heat  conductivity,  emissivity,  albedo  and  density 

Used  in  calculation  of  DEPTH  when  converting  new  snow 
to  old  snow,  cm. 
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AB 

ACL(8) 

ALPH(IX) 

AO 

APRM 

ATERM 

AVABSOR 

AVATERM 

AVDTERM 

AVGBR 

AVHTERM 

AVSOL 

AVSTR 

B 

BBB(J,I) 

BC 

BCL(8) 

BEP 

BK 

BPRM 


Table  II  :  Variable  definitions  :  TSTM 

MUGGE-MOLLER  ABSORPTION  FUNCTION. 

COEFFICIENT, a,  DEPENDANT  ON  CLOUD  TYPE.  USED  IN  THE 
CALCULATION  OF  CTF. 

THERMAL  DIFFUSIVITY  OF  LAYER  IX  IN  CM**2/MIN 
ATMOSPHERIC  ALBEDO  FOR  RAILEIGH  SCATTERING. 
FACTE*TEMP**3  IN  CAL/CM**2-MIN-C 

ENERGY  CONTRIBUTED  BY  ATMOSPHERIC  IR  EMISSION 
CAL  CM**2-MIN 

AVERAGE  DAILY  SHORTWAVE  RADIATION  ABSORBED  AT  THE 
SURFACE.  Win-2. 

AVERAGE  DAILY  INCOMING  LONGWAVE  RADIATION , Wm- 2 . 

AVERAGE  DAILY  EVAPORATIVE  HEAT  FLUX,  Wm-2. 

AVERAGE  DAILY  GREYBODY  RADIATION,  ie.  REFLECTED 
LONGWAVE  RADIATION,  Wn»-2. 

AVERAGE  SENSIBLE  HEAT  FLUX.  Wm-2. 

AVERAGE  DAILY  INCOMING  SHORTWAVE  RADIATION,  Wm-2. 

AVERAGE  SNOW  TEMPERATURE,  oC. 

HEAT  CONDUCTIVITY  OF  SURFACE  CAL/CM**2-MIN-C 

Y  INTERCEPT  OF  LINEAR  EQUATION,  USED 
FOR  TABLE  INTERPOLATION. 

CONSTANT  USED  IN  THE  CALCULATION  OF  WATER  VAPOUR 
PRESSURE. 

COEFFICIENT,  b,  DEPENDANT  ON  CLOUD  TYPE.  USED  IN  THE 
CALCULATION  OF  CTF. 

BOTTOM  BOUNDARY  THERMAL  IR  EMISSIVITY. 

BOTTOM  SURFACE  GEOMETRIC  SHAPE  IN  FRACTION(0 . 0- 1 . 0) 
HEAT  CONDUCTIVITY  OF  BOTTOM  BOUNDARY  LAYER 


/ cont. . . 
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BTERM  ENERGY  CONTRIBUTED  BY  INSOLATION  AFTER  ADJUSTMENT  USING 

SURFACE  ABSORPTIVITY.  IN  CAL/CM**2-MIN 

CC  CLOUD*CLOUD . 

CCOEF  PART  OF  CALCULATION  OF  INTERFACE  VALUES. 

CLR(8)  COEFFICIENT  DEPENDANT  ON  CLOUD  TYPE,  USED  IN  THE 
CALCULATION  OF  INCOMING  LONGWAVE  RADIATION. 

COEI  USED  IN  CALCULATION  OF  SCF,  DETERMINED  BY  THE  RICHARDSON 

NUMBER. 

COE2  USED  IN  CALCULATION  OF  SCF.  DETERMINED  BY  THE  RICHARDSON 

NUMBER. 

CLOUD  CLOUD  COVER  IN  FRACTION  OF  0.1 -1.0. 

CP  SPECIFIC  HEAT  OF  DRY  AIR  AT  CONSTANT  PRESSURE. 

CTEMA  AIR  TEPERATURE  IN  DEGREES  KELVIN,  USED  IN  CALCULATION  OF 

EVAPORATIVE  HEAT  FLUX. 

CTF  CLOUD  ADJUSTMENT  FACTOR. 

DAIMX  MATRIX  HOLDING  TSTM  OUTPUT  VALUES  FOR  MANIPULATION, 

day  JULIAN  DAY  USED  IN  SOLVING  INSOLATION 

DCOEF  PART  OF  CALCULATION-OF-INTERFACE  VALUES. 

DDDT  dT/dTg.  PART  OF  UPPER-BOUNDARY-CALCULATION. 

DECL  SOLAR  DECLINATION  ANGLE 

DELT  TIME  STEP  IN  HOURS 

DEPTH (lY)  MATRIX  HOLDING  LAYER  DEPTHS. 

DEPTH(450)  DEPTH  OF  SNOW  OR  SOIL  PROFILE. 

DIST  DEPTH  IN  CM  OF  INITIAL  SOIL  PROFILE  AT  WHICH 

CORRESPONDING  SOIL  TEMPERATURE  IN  DEGREE  C  IS 
INTERPOLATED.  (TABLE  5). 

DNEW  D  VALUE  USING  LATEST  ESTIMATE  OF  GROUND  TEMPERATURE. 

DOWNIR  INCOMING  LONGWAVE  RADIATION,  CM-2  MIN-1. 

DPRM  HEAT  FLUX  IN  CAL/CM**2-MIN  AT  BOTTOM  BOUDARY  OR 


/cont. .. 
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TEMPERATURE  IN  RANKINS  AT  BOTTOM  BOUNDARY. 


DPRMO 

DPRMl 

DTERM 

DUST 

ELF 

EPSN 

ES 

ESAT(T) 

E(T) 

EX 

r2 

FACTA 

FACTD 

FACTE 

FACTH 

FK(IX) 

FMM(J,I) 

G 

GTERM 


TEMPERATURE  OF  BOTTOM  MATERIAL  IN 
DEGREE  CELSIUS. US ED  WHEN  LFLUXY-0 

HEAT  FLUX  OF  BENEATH  BOTTOM  MATERIAL. 

IN  CAL/CM**2-MIN,  USED  WHEN  LFLUXY  NOT 
EQUAL  0. 

ENERGY  LOSS  DUE  TO  EVAPORATION. 

ATMOSPERIC  DUST  IN  POUNDS/CUBIC  CENTIMETERS 
(LBS/CC)USED  IN  SOLVING  INSOLATION. 


LATITUDE  IN  RADIANS 
EMISSIVITY  OF  SURFACE  MATERIAL 

SATURATED  VAPOUR  PRESSURE  ^iT  TEMPERATURE  T. 

VAPOUR  PRESSURE  AT  TEMPERATURE  T. 

USED  IN  CALCULATION  OF  SCF.  DETERMINED  BY  THE  RICHARDSON 
NUMBER. 


PART  OF  CALCULATION -OF- BOTTOM  BOUNDARY  VALUES. 
SIGMA*EPSN 

FACTD-SIGMA*BK*BEP*TR**4  USED  IN  BOTTOM  BOUNDARY 
CALCULATION  WHEN  THERE  IS  AIRSPACE  BENEATH  THE  BOTTOM 

FACTE-SIGMA*BK*BEP  USED  IN  BOTTOM  BOUNDARY  CALCULATION 
WHEN  THERE  IS  AIRSPACE  BENEATH  THE  BOTTOM 

USED  IN  SOLVING  CONVECTION  TERM  (HTERM) 

( 1000 . 0/PRESS ) **0.286 

HEAT  CONDUCTIVITY  OF  LAYER  IX  IN  CAL/MIN-CM-K 
SLOPE  OF  LINEAR  EQUATION, USED  FOR  TABLE  INTERPOLATION 

ACCELERATION  DUE  TO  GRAVITY. 

GROND  HEAT  FLUX,  Wm-2 


/cont. . . 
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HEADER  72  CHARCTER  INPUT  VARIABLE  USED  TO  PRINT 
COMMENTS  ON  OUTPUT. 

HTERM  ENERGY  LOSS  OF  GAIN  DUE  TO  CONVECTION  CAL/CM**2 - MIN 

lABSOR  INCOMING  SHORTWAVE  RADIATION,  Wm-2. 

lATERM  INCOMING  LONGWAVE  RADIATION.  Win-2. 

ICNT  TSTM  INTERNAL  ITERATION  COUNT. 

IDTERM  EVAPORATIVE  HEAT  FLUX,  Wni-2. 

lEFSWT  SWITCH  WHEN  -0  WILL  PRINT  OUTPUT  ONLY  AT  SPECIFIED  TIME 
IF  NOT  -0  WILL  PRINT  OUTPUT  AT  EVERY  ITERATIONS. 

lEOF  SET  FROM  0  TO  1  WHEN  AN  EOF  IS  ENCOUNTERED.  USED  TO 

TERMINATE  PROGRAM 

IHTERM  SENSIBLE  HEAT  FLUX,  Wni-2. 

Ill  COUNT  OF  NUMBER  OF  INTERNAL  ITERATIONS  IN  CALCULATION 

OF  UPPER-BOUNDARY-VALUES. 

IMATL  BACKWARD  COUNTER  OF  LAYERS.  STARTING  WITH  THE  NUMBER 

OF  LAYERS. 

INTR(IX)  BEGINNING  SUB-LAYER  DEPTH  NUMBER  FOR  LAYER  NUMBER  IX 

IPRNT  BACKWARD  COUNTER  SET-NPRNT.  WHEN  EQUAL  TO  1  OUTPUT  IS 

PRINTED . 

ITER  ITERATION  COUNTER  USED  IN  FINITE  DIFFERENCE  CALCULATION 

OF  HEAT  FLOW  EQUATION 

ITIME  BACKWARD  COUNTER  INITIALIZE  AS  TOTAL  TIME  STEPS  IN  HOUR 

IX  LAYER  NUMBER  STARTING  WITH  TOP  LAYER 

lY  SUB -LAYER  DEPTH  NUMBER 

JMAX  THE  TOTAL  NUMBER  OF  SUB -LAYERS 

KIO  TSTM  INTERNAL  ITERATION  COUNT. 

KSQ  VON  KARMAN'S  CONSTANT  SQUARED. 

KTEMPA  AIR  TEMPERATURE  IN  DEGREES  KELVIN. 

KTEMPG  SURFACE  TEMPERATURE  IN  DEGREES  KELVIN. 

LAYERS  COUNT  OF  THE  PROFILE  LAYERS. 

LAT  LATITUDE  USED  IN  SOLVING  INSOLATION 


/oont. .. 
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LFLUXY 


LN 

M 

MAX(J) 

NCLOUD 

NIT 

NIPTS 

NOMATL 

NPRNT 

NTABL 

NX(IX) 

OUTCD 

PI 

PPTN 

PRESS 

PTOIE 

QSAT(T) 

Q(T) 

REP 

RK 

RHOA 


INPUT  BOTTOM  BOUNDARY  DATA  CONTROL  SWITCH.  IF-0,THERE 
IS  NO  HEAT  FLUX  THROUGH  BOTTOM  OF  MATERIAL, IF  NEGATIVE 
THERE  IS  NO  AIR  SPACE  BENEATH  BOTTOM  MATERIAL, IF  POSIT¬ 
IVE  THERE  IS  AIR  SPACE  BENEATH  BOTTOM  MATERIAL. 

DUMMY  VARIABLE  TO  READ  LINE  NUMBER  FROM  INPUT  FILE 


SECANT  OF  SOLAR  LENITH  AliGLE  IN  RADIANS 

NUMBER  OF  INPUT  TABLE  VALUES  USED  IN  TABLE 
INTERPOLATION  MODULE. 

CLOUD  TYPE  INDEX  NUMBER  (1-9)  USED  IN 
SOLVING  INSOLATION, INFRARED  EMISSION. 

NUMBER  OF  ITERATIONS  USED  IN  HEAT  FLOW  CALCULATIONS 

NUMBER  OF  COMPUTATIONAL  POINTS  WITHIN  PROFILE. 

NUMBER  OF  MATERIAL  LAYERS  USED  IN  SOLVING  HEAT  FLOW 

NUMBER  OF  TIMES  OUTPUT  TIME  PRINT  FREQUENCY  IS  DIVISIBL 
BY  TIME  STEPS.  USED  TO  DETERMINED  WHEN  TO  PRINT  OUTPUT. 

TABLE  NUMBER 

NUMBER  OF  SUBLAYER  OF  EACH  LAYER, NX(IX)-THK(IX)/SFRQ(IX 

OUTPUT  MANIPULATION  VARIABLE. 


3.141593 

PRECIPITATION,  MM  OF  WATER. 

ATMOSPHERIC  PRESSURE  IN  HILLIBAR(HB) 

USED  IN  SOLVING  INSOLATION 

BEGINNING  TIME  OF  OUTPUT-TOTAL  NUMBER  OF  HOURS  MINUS  24 
USED  IN  PRINT-OUTPUT  MODULE  2 


SATURATED  SPECIFIC  HUMIDITY  AT  TEMPERATURE  T. 
SPECIFIC  HUMIDITY  AT  TEMPERATURE  T. 


EMISSIVITY  BENEATH  AIRSPACE 
RELATIVE  HUMIDITY,  % 

AIR  DENSITY,  CALCULATED  AS  PART  OF  SENSIBLE  HEAT 
CALCULATION. 


/ cont. . , 
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RHOC(IX)  FK(IX)/ALPH(IX)  IN  CAL/CM**2-K 


RI  RICHARDSON  INDEX  NUMBER  USE  IN  SOLVING  CONVECTION 

ENERGY  LOSS. 

RK  SURFACE  BENEATH  AIRSPACE  GEOMETRIC 

SHAPE  IN  FRACTION  (0.0  -  1.0) 

RR(IX)  RR(IX;-i)ELT/SFREQ**2.  (PART  OF  HEAT  FLOW  EQUATION) 

SAZ  SOLAR  AZIMUTH  IN  RADIANS.  SAZ-ATAN( -COS(DECL)*SIN(TIMER 

( COS ( ELF*SIN ( DECL) - S IN ( ELF) *COS (TIMER)  ) ) 

SCF  USED  IN  CALCULATION  OF  SENSIBLE  HEAT  FLUX,  DETERMINED  BY 

THE  RICHABDSN  NUMBER. 

SFRQ(IX)  VERTICAL  GRID  SPACING  IN  CM  IN  EACH  LAYER  IX  IN  CM**2/M 

SICF  INSOLATION  ADJUSTMENT  DUE  TO  ZENITH  ANGLE, SURFACE  SLOPE 

AND  SURFACE  ASPECT  ANGLE.  SICF-COS(Z)*COS(SLOPE)+SIN(Z) 

S  IN  ( SLOPE ) *COS ( SAZ - SURFAC ) 

SIGMA  STEPHAN- BOLTZMAN  CONSTANT,  a.l2E-ll 

SLOPE  SURFACE  SLOPE  IN  DEGREES  WITH  HORIZONTAL-0  DEGREE, 

USED  IN  SOLVING  INSOLATION 

SMALLA  ABSORBTIVITY  OF  SURFACE  MATERIAL 

SPEED  WIND  SPEED  IN  CM/SEC 

ST0R(1,IY)  ESTIMATE  SUB-LAYER  TEMPERATURE  IN  DEGREE  RANKINE 

STOR(2,IY)  FK;HEAT  CONDUCTIVITY  OF  SUB-LAYER  lY  IN  CAL/MIN-CM-K 

STOR(3,IY)  RHOC,FK/ALPH  IN  CAL/CM**2.k 

ST0R(4,IY)  CONSTANT  DIMENSIONLESS. 

STOR(5,IY)  INITIAL  SOIL  TEMPERATURE  IN  DEGREE  RANKINS 
OF  INITIAL  SOIL  PROFILE 

STOR(6,IY)  SAME  AS  STOR(2,IY) 

STOR(7,IY)  SAME  AS  STOR(3,IY) 

SUN  CALCULATED  INSOLATION  VALUE. 

SURFAC  SURFACE  AZIMUTH  IN  DEGREE  WITH  SOUTH  -0  DEGREE, US ED 
IN  SOLVING  INSOLATION 

SAME  AS  TIME 

AIR  TEMPERATURE  IN  DEGREE  RANKINE 
TOTAL  DAILY  INCOMING  LONGWAVE  RADIATION,  Woi-2. 


T 

TA 

;  TATERM 

r 

I 


IV-7 
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TABSOR 

TAG 

TAK 

TAL 

TB 

TD 

TDTERM 

TEML 

TEMR 

TFRQ 

TGBR 

THK(IX) 

THTEEM 

TIME 

TIMER 

TOTTIM 

TPRNT 

TR 

TSK 

TSOL 

TYME 

WATER 

WET 


TOTAL  DAILY  SHORTWAVE  RADIATION  ABSORBED  AT  THE 
SURFAGE,  Wni-2. 

AIR  TEMPERATURE  IN  DEGREE  GELSIUS 

air  TEMPERATURE  IN  DEGREE  KELVIN 

GONSTANT  USED  IN  THE  GALCULATION  OF  WATER. 

THERMAL  GONDUGTIVITY  OF  BOTTOM  MATERIAL 
GAL/CM**2-DEG  G-MIN 

DEW  POINT  TEMPERATURE  USED  IN  THE  GALCULATION  OF  WATER. 
TOTAL  DAILY  EVAPORATIVE  HEAT  FLUX,  Wni-2. 

SURFACE  TEMPERXTURE  OF  MATERIAL  IN  DEGREES  RANKINE. 
BOTTOM  LAYER  TEMPERATURE  OF  MATERIAL  IN  DEGREES  RANKINE. 
TIME  STEP  IN  MINUTES  USED  IN  SOLVING  HEAT  FLOW 
TOTAL  DAILY  REFLECTED  RADIATION,  Win.2. 

LAYER  THICKNESS  IN  CM  OF  LAYER  IX 
TOTAL  DAILY  SENSIBLE  HEAT  FLUX,  Wm.2. 

TIME  IN  HOURS  IN  WHICH  MATERIAL  TEMPERATURES 
ARE  ESTIMATED 

SUN'S  HOUR  ANGLE  IN  RADIANS 

TOTAL  NUMBER  OF  24  HOUR  REPETITIONS  USED  IN  SOLVING 
HEAT  FLOW 

OUTPUT  TIME  PRINT  FREQUENCY  IN  MINUTES 
TEMPERATURE  OF  AIRSPACE  BENEATH  BOTTOM  MATERIAL. 

MATERIAL  SUB -LAYER  TEMPERATURE  IN  DEGREES  KELVIN 
TOTAL  DAILY  INCOMING  SHORTWAVE  RADIATION,  Wm-2. 

TIME  IN  HOURS  USE  INSOLATION  CALCULATION 

THE  AMOUNT  OF  PRECIPIPAL  WATER  IN  MILLIMETERS 
(MM)  USED  IN  SOLVING  INSOLATION. 

MOISTURE  CONTENT  OF  SURFACE  MATERIAL 
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XXX(J,1) 

XXX(J,2) 

XXX(J,6) 

XXX(J.3) 

XXX(J,4) 

xjacj.s) 

YY?(J.l) 

YYY(J,2) 

YYY(J,6) 

YYY(J,3) 

YVY(J.4) 

Z 

ZA 

ZASH 

Z2A. 

ZZB 


TIME  IN  HOURS  FOR  TABLE  1  (AIR  TEMPERATURE) 

TIME  IN  HOURS  FOR  TABLE  2  (RELATIVE  HUMIDITY) 

TIME  IN  HOURS  FOR  TABLE  6  (WIND  SPEED) 

TIME  IN  HOURS  FOR  TABLE  3  (AMOUNT  OF  CLOUD  COVER) 

TIME  IN  HOURS  FOR  TABLE  4 

DEPTH  IN  CENTIMERERS  FOR  TABLE  (5) (TEMPERATURE  PROFILE) 
INITIAL  TEMPERATURE  PROFILE  IN  DEGREES  CELSIUS 


AIR  temperature  IN  DEGREES  CELSIUS  TABLE  1 

RELATIVE  HUMIDITY  INPUT  FRACTION, USED  IN  TABLE  2  SOLVIN 
INFARED  EMISSIONS, (ATERM)  EVAPORATIVE  HEAT  LOSS  (DTERM) 

WIND  SPEED  INPUT  IN  METERS/SECOND  AND  CONVERTED  TO 
CENTIMETER/SECOND  TABLE  6 

AMOUNT  OF  CLOUD  COVER  IN  FRACTION  (0  TO  1)  USED  IN 
SOLVING  INSOLATION  TABLE  3  INFARED  EMISSION  (ATERM) 

INSOLATION  IN  CAL/CM**2-MIN, IF  0.0  AT  1200  HOURS, 
INSOLATION  VALUES  WILL  CACULATED  TABLE  4 

SOLAR  ZENITH  ANGLE.  Z-SIN(DECL)*SIN(ELF)+COS(DECL)* 

COS (ELF)*C0S (TIMER) 

SHELTER  HEIGHT  IN  CENTIMETERS (CM) 

SHELTER  HEIGHT  IN  CENTIMETRES (CM) , 

SURFACE  TEMPERATURE  OF  MATERIAL  IN  DEGREE  RANKINE 

BOTTOM  LAYER  TEMPERATURE  OF  MATERIAL  IN  DEGREE  RANKINE 
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V-1 


K\tr 


c  ************  program  snomo.f  ************* 

c 

c  **SNOMO.F  IS  A  DISTRIBUTED  SNOWMELT  MODEL.  THE  MODEL  OPERATES  IN 
c  CELLS,  EACH  CELL  IS  A  SUBDIVISION  OF  THE  CATCHMENT  THAT  IS  BEING 
c  MODELLED  AND  IS  ASSUMED  TO  HAVE  A  HOMOGENEOUS  VEGETATION  COVER, 
c  ALTITUDE  AND  ASPECT.  THE  SNOWMELT  RATE  IS  DERIVED  FROM  THE 
c  CALCULATION  OF  THE  SNOWPACK  ENERGY  BUDGET.  THE  EFFECTS  OF  MELT 
c  AND  SNOWFALL  ON  THE  SNOWPACK  ARE  MODELLED, 
c 
c 

c  **Initlalise  variables,  constants,  matrices  etc.. 

c  **Data  is  input  from  2  files.  Veg.dat  contains  only  vegetation 

data 

c  and  is  called  when  iveg  is  1.  Snomo.dat  contains  daily 
meteorological 

c  data, physical  constants  and  computational  controls  for  both 
vegetation 

c  and  non-vegetation  situations, 

c. 

common/matrixl/xxx(30, 10) ,yyy(30, 10) 
cominon/matrix2/depmxl(2 , 2)  ,depmx2(3 , 2) 
coninion/maCrix3/thkmxl(l,6) ,  thla]ix2(2, 6) 
common/tstml/press , ncloud,za, slope , surfac , day, lat ,wec 
common/ tstm2/nomaCl 

common/ ts tm3/Cottim, tfrq, tprnt,nipts 
common/vegi/iveg, sigf , state , epf , fola ,hf ol 
common/rain/airtp 
common/s  frq/s  frq ( 6 ) 
common/phycons/lheatf , hcapsn , hcapw 
real  area , sur f tp , snpptn , snvol , critdp , 

&  enavm, del ts , hcapsn, densn, rstpc, denw, hcapw, Iheatf , 

&  avsndp, depth, coefb,pcerm,netrad, 

&  avbottp ,pptn, dayl , lat , tgbr . tsol , tabsor , taterm , 

&  thterm, tdterm, vsoil(l ,6) ,vnewsn(l , 6) , voldsn(l , 6 ) , 

&  wel , depthl ,mrwe ,mrsn,mrlwe ,mrlsn, enavml 
integer  i,j , n,kl ,k2 ,k3 ,klO ,kontrol 

data  denw, hcapw, hcapsn, Iheatf/lOOO .0,4.21,2.09,0. 334e+6/ 

data  tottim, tfrq , tprnt/1 , 1 . ,60./ 

nin-10 

nveg-11 

nout-9 

open(unit-nin, f lie-' snomo .dat' , status-' old' ) 
rewind  nin 

open(unit-nveg, file-' veg. dat' , status-' old' ) 
rewind  nveg 

open(unlt-nout , f ile-' SNOMO. RES' ) 
rewind  nout 


read(10 ,*) ( (vnewsn(i, J ) . j-1 , 6) , i-1 , 1) 
read(10,*)((voldsn(i.j) ,J-1.6) .1-1,1) 
read(10,*)((vsoil(l,j) ,j-l,6) .1-1,1) 
read(10 ,*) press .ncloud.za 
read(10 , *) slope , surfac , dayl , lat 
read(10,*)lveg 
If (Iveg. eq. 1) then 

read( 11 , *) s Igf , s  Cate , epf . f ola , hf o 1 
end  If 

read(10  ,*)_area,  crlcdp ,  deles  .coefb  .wet 
read(10,*)solldp,sollcp,ssolltp 
read(10,*)n 
c 

c**  n  number  of  days  to  be  modelled  le .  no.  of  days  in  data  file 

c  minus  one . 

c**  xxx(l,8)-  Julain  date 

c  xxx(l.l)-  observation  time,  24hr  clock 

c  yyy(l,l)-  air  temperature,  oC 

c  yyy(l,2)-  relative  humidity,  % 

c  yyy(l,3)-  cloud  cover,  0-1 

c  yyy(l,6)”  wind  speed,  m/s 

c  yyy(1.7)-  precipitation,  mm  water 

c 

read(10,*)xxx(l, 8) ,xxx(l , 1) ,yyy(l , 1) ,yyy(l. 9) ,yyy(l .  2)  , 

&  yyy(l,3), 

&  yyy(i.6) .yyy(i.7) 

read(10,*)xxx(2 ,8) ,xxx(2 , 1) .yyy(2 , i) ,yyy(2 , 9) ,yyy (2 , 2)  , 

&  yyy(2,3), 

&  yyy(2,6) ,yyy(2,7) 

write(9 ,*)' input  matrix  ' 

write(9 ,*)xxx(l,8) ,xxx(l, 1) ,yyy(l , 1) .yyy(l , 9) ,yyy(l . 2) , 

&  yyy(l,3), 

&  yyy<1.6) ,yyy(l,7) 

wrice(9,*)xxx(2,8) .xxx(2.1) ,yyy(2,l) .yyy(2,9) ,yyy(2,2)  . 

&  yyy(2,3), 

&  yyy(2,6) ,yyy(2.7) 

c 
c 

do  20  1-1,2 
do  10  j-1.6 
thkmx2(i,j)-0.0 
depmx2(l, j )-0.0 
10  continue 

20  continue 

konCrol-0 
snvo 1-0.0 
depth-0 . 0 
kl-0 
k2-0 


k 


c 

c  **The  dally  computational  loop  is  do  5001 
c 

do  5001  1-1, n 
kl-kl+1 

kelvlii-yyy(l ,  l)+273 . 15 
pptn-yyy(l,7) 
airt  -yyy(l,l) 
ddate-jcxx(l,  8) 
time2-xxx(2 , 1) 
c 

c**If  statement  to  determine  if  is  the  first  day  of  model  run 
c  and  if  the  has  been  any  snowfall, 
c 

if (kl . eq. 1) then 
write(_9.*) 'FLAG  20' 

writs(?,*)'l  PPTN.AIRTP-  ' , snpptn, airtp 
k2-0 

if (pptn. gt . 0 . 0 . and. airtp . It .0 . 0) then 
goto  50 
end  if 
end  if 
c 

c**If  statement  to  determine  if  there  was  a  snowcover  present  in  the 
c  previous  iteration, 
c 

if(k3.eq.l)then 

k2-0 

goto  57 

else 

depth-0 . 0 
goto  50 
end  if 
c 

50  write (9,*) 'FLAG  21' 

c**Calculate  snow  precipitation. 

snpptn-(pptn/vnewsn(l , 6) )/10 
write(9 ,*) ' snpptn (cm) .pptn(mm)-  ' , snpptn, pptn 
dep  th-snpp  tn+dep  th 
c 

c**Snow  present:  set  depth  matrix  using  surface  soil  and  air 
c  temperatures  as  snow  surface  and  basal  temperatures, 

c  snow  state  is  'new' . 

c**Snow  absent:  set  depth  matrix  using  the  known  soil  depth,  soil 
temperature  at  depth  and  air  temperature  as  the 
soil  surface  temperature. 


c 

c 

c 


If (depth. ge . critdp)chen 
write (9,*) 'FLAG  22' 
nomat 1-1 
nlpts-2 

depmxlCl , l)-0 . 0 
depmxl(l , 2)-airtp 
depmxl(2 , l)-depth 
depmxl(2 , 2)-ssoiltp 
do  52  1-1,2 

write  (9, 555)  (depmxK  i  ,  j  )  ,  j-1 , 2) 

52  continue 

do  54  j  -1 , 6 

thkmxl(l, j )-vnewsn(l, j ) 

54  continue 

thkmxl(l , 1) -depth 

write (9.*) 'THICKNESS  MATRIX;' 

write (0  ,  *)  (  <  thkmxKi ,  j  )  ,  j-1 , 6)  ,  i-1 , 1) 

write(9,*) 'FLAG  23' 

k25-l 

k3-l 

k2-l 

goto  57 

else 

wrlte(9,*)'FLAG  24' 
nomatl-l 
nipts-2 
do  56  j-1 , 6 
rhkmxi(l , j )-vsoil(l , j ) 

56  continue 
thknixl(l ,  l)-soildp 
depmxl( 1 , l)-0 . 0 
depmxl (1 , 2)-airtp 
depmxl (2 , 1) -soil dp 
depnixl(2 , 2)-solltp 
k3-0 

k2-0 
end  if 
c 

c**Calculatlon  of  sfrq  using  subroutine  csfrq. 
c  sfrq(ix)  is  the  vertical  grid  spacing  in  cm  in  each  layer  ix 
c  in  cm**2/m.  sfrq  is  variable  because  layer  size  is  variable, 
c 

57  call  csfrq(l,n) 
c 

c**2-layered  pack:  set  depth  matrix  for  1-layered  pack  to  0.0 
c  1-layered  pack:  set  depth  matrix  for  2-layered  pack  to  0.0 
c 


if (nomatl . eq . 1 ) then 
do  59  i-1,3 
do  58  j-1.2 
depmx2 ( 1 , j ) -0 . 0 

58  continue 

59  continue 
else 

do  61  i-1.3 
do  60  j-1.2 
depmxl ( i , j ) -0 . 0 

60  continue 

61  continue 
end  if 

c 

do  62  i-1,3 

write(9 , 555) (depmx2(i , j ) , j-1 , 2) 

62  continue 

do  63  i-1,2 

write (9 ,555) ( depmxl (i.j ) , j-1 , 2) 

63  continue 
c 

c**Set  some  variables  to  zero  before  calling  TSTM 
c  The  write  statements  are  a  program  developement  check, 
c  The  format  statements  are  for  the  output. 
wrlte(9,*)'FLAG  25' 
write(9 , *) ' sfrq-  ',sfrq 
klO-0 
avstp-0.0 
avsol-0 . 0 
avgbr-0 . 0 
avabsor-0 . 0 
ava term-0.0 
avhterm-0 . 0 
avdterm-0 . 0 
avbottp-0.0 

wrlte(9.*) 'A1_ED0  1-  ' . thkmxl(l , 5) 
wrice(9,*) 'ALBEDO  2-  ' . thkmx2(2 . 5) 
write(9,*) '**DATA  MATRIX  PASSED  TO  TSTM**' 
write (9 , 64)xxx(l, 8) ,xxx(l, 1) .yyy(l , 1) ,yyy(l , 2) ,yyy(l , 3) , 
&  yyy(1.6) .yyy(l,7) 

wrlte(9,64)xxx(2 , 8)  ,xxx(2, 1)  ,yyy(2 , 1)  ,yyy(2 , 2)  ,yyy(.2 , 3) , 
&  yyy(2,6) ,yyy(2,7) 


64  format(7f6.1) 

310  FORMAT (1 IX. 'TOTAL  GRAYBODY  EFFECTIVE  GROUND  FOLIAGE’, 

6.  4X, 'SOLAR') 

320  FORMAT (14X. 'RADIANCE  TEMP' . lOX, 'TEMP  TEMP’ 

&  ,4X, 'INSOLATION') 

330  F0RMAT(5X. 'HR' ,7X, ' (W/M**2)' ,9X, ' (C)' .IIX, ' (C) ’ 

&  ,6X,'(C)  (W/M**2)') 

340  F0RMAT(9X,' . REFL-NREFL - REFL - NREFL’  ,30(1H-)) 

270 

FORMAT(3X,F5.2,5X,I4,2X.I4,3X,F6.1.2X.F6.1,3X,F6.1.3X,F6.1,7X,I4) 

355  FORMAT ( 5 7X, 'SENSIBLE  LATENT') 

360  FORMAT (5X, 'HRS  SURFACE  GRAYBODY  SOLAR  SURFACE  ATMOS’, 

&  ’  IR  HEAT  HEAT') 

370  FORMAT (1 IX. 'TEMP  RADIANCE  INSOLATION  ABSORP  EMISSION’, 

&  '  LOSS  LOSS’) 

380  FORMAT (1 IX, 'DEG  C' ,23(1H- ) . ' (W/M**2) ' . 24(1H- ) ) 


390  FORMAT(3X. 'TIME' .18X, 'DEPTH.  TEMPERATURE ' /4X HR 2 IX CM ' , 

&  9X, 'C'/2X,65(1H-)) 
if ( iveg . eq . 0) then 
WRITE(9.355) 

WRITE(9,360) 

WRITE(9,370) 

WRITE(9,380) 

else 

WRITE(9,310) 

WRITE(9,320) 

WRITE(9,330) 

WRITE(9,340) 

endlf 

write(6,*) '**ENTERING  TSTM  SUBROUTINE**' 
c 

c**Call  the  subroutine  TSTM.  TSTM  calculates  the  shortwave,  longwave, 
c  latent  heat  and  sensible  heat  energy  exchanges  at  the  surface  of 
the 

c  snowpack  and  the  heat  transport  through  the  pack.  TSTM  calculates 
most 

c  of  the  variables  used  in  the  calculation  of  snowmelt. 

c  In  SNOMO  TSTM  can  model  either  snow  or  non- snow  covered  ground. 

c 

call  tstm( i , ' Y' , ' N' ,kl0 , avstp .avsol , avgbr , avabsor , avaterm , 

&  avhterm, avdterm, k2 , avbottp , tgbr , tsol , tabsor , 

&  taterm, thterm, tdterm) 

write(6,*) '**EXITED  TSTM  SUBROUTINE**' 
write(9,*) 'avstp, avgbr, avsol  etc.  ' 

wrlte(9 ,*)avstp , avgbr , avsol, avabsor , avaterm, avhterm, avdterm 
write (9 ,*)' average  base  temp.-  '.avbottp 
c 

if  (k3 . eq . 0) then 

write(9 ,*) ' TSTM  with  SOIL  characteristics  was  called' 

goto  9000 

else 

k3-l 

wrlte(9 , *) ' TSTM  with  SNOW  characteristics  was  called' 
goto  65 
end  if 
c 

65  pterm-0.0 

c 

c**Call  subroutine  cmelt.  This  calculates  meltrate  and  volume  of  snow 
c  melted.  Clmelt  is  identical  except  that  dally  averages  of  the 
c  energy  exchanges  are  used,  not  totals.  Meltrate  is  calculated  in 
c  cm  of  snow  and  mm  of  water  equivalent, 
c 

call  cmelt(l,n,enavm,netrad,pterm, tsol, tgbr  tabsor, 

&  taterm, thterm, tdterm, pptn,densn,mrwe ,mrsn, gterm, denw, depth, 

&  avstp , avbottp) 

wrlte(9,*)'l  ENAVM-  ' ,enavm 

wrlte(9,*)'l  HRWE(mm. water) ,MRSN(cm. snow)-  '.mrwe.mrsn 
wrlte(9, *)’**' 

call  clmelt( 1 ,n, enavml .netrad, pterm, avsol , avgbr , avabsor , 

&  avaterm,avhterm,avdterm,ppcn,densn,mrlwe ,mrlsn, gterm, denw, 

&  depth, avstp, avbottp) 

write (9,*) '2  ENAVMl-  '.enavml 

wrlte(9,*)'2  MRlWE(mm. water) ,MRlSN(cm. snow)-  ' ,mrlwe .mrlsn 


write(9,*)'  ' 


c**Rain  and  melcing  pack:  call  subroutine  raininp.  Melcrace  has  to  be 
c  recalculated. 

c  Rain  and  frozen  pack:  call  subroutine  ralnfp 
c 

if (pptn. gt. 0. 0. and. airtp. gt. 0. 0) then 
if (mrwe . gt . 0 . 0) then 
call  ralnnp(i,n,ppcn,pcena,denw) 

call  cmelc(l,n,enavni,necrad,pcenn,  tsol,  tgbr,  cabsor, 

&  taterm, thterm, tdtenn, pptn, densn, mrwe ,mrsn, gcerm, denw, 

&  depth, avstp.avbottp) 

write(9,*)'8  ENAVM-  '.enavm 

write(9,*)'8  MRWE(nBn.  water)  .MRSM(cni.  snow)-  ',mrwe,mrsn 
wrlte(9 ,*) '**' 

call  clmeltCl , n.enavml .netrad, pterm, avsol , avgbr , avabsor , 

&  avaterm.avhterm, avdterm.pptn, densn.mrlwe ,mrlsn, gterm, 

&  denw.depth.avstp.avbottp) 

write(9,*)'9  ENAVMl-  ' .enavml 

wrlte(9,*)'9  MRlWE(nmi.  water)  ,MRlSN(cm.  snow)-  '  ,mrlwe,mrlsn 
write(9,*)'  ' 
goto  66 
else 

write(9.*)'**RAIN-0N- FROZEN  PACK  CALLED**' 
call  ralnfp ( i , n , pp tn , r s  tpc , denw , avs  tp , dep  th , densn) 
end  if 
end  if 
c 

c**Calculate  snow  precipitation  if  is  not  the  first  SNOMO  run  or 
c  the  first  day  with  snowfall. 

66  if(k2.eq.0)chen 

if (pptn . gt . 0 . 0 . and. airtp . It . 0 . 0) then 
snpptn-(ppcn/vnewsn(l , 6) )/10 
wrice(9,*)'2  pptn(mni)-  '  ,pptn 
wrlte(9,*)'2  snpptn(cm)-  ' .snpptn 
else 

snpptn-0 . 0 

wrlte(9,*)'3  snpptn(cni)-  '.snpptn 
end  if 
goto  68 
end  if 
c 

c**Manlpulatlon  of  the  snowpack  to  allow  for  the  effects  of  snowfall, 
c 

68  if(k25.ge.3)then 

if(snpptn.gt.0.0)then 
if (nooatl . eq . l)chen 
write(9,*) ’chkntxl(l.l)-  ' ,  thkmxKl,  1) 

wo 1- ( thknxl ( 1 , 1 ) * 10 ) *vnewsn ( 1 , 6 ) 
depthl-(wel/voldsn(l . 6) )/10 
depth-depchl 

wrlto(9 ,*) ' 10  wel-  ’ pWol 
wrlte(9 ,*) ' 10  depth-  '.depth 
write(9,*) 'FLAG  2' 
goto  125 
end  if 


write (9.*) 'FLAG  3’ 
c  depth-depth-snpptn 

wel-(chknix2(l ,  l)'*10)*vnewsn(l ,  6) 

depthl-(wel/voldsn(l . 6) )/10 
depth-depchl 
depth-Chknix2  ( 2 ,  l)+depth 
125  nomatl-2 

nipcs-3 
do  150  J-1,6 
thknix2 ( 1 ,  j  )-vnewsn(l,  j  ) 
thlciiix2 ( 2 ,  j  )-voldsn(l,  j ) 

150  continue 

thknix2  (1,1)  -snpp  tn 
thknix2  (2,1)  -dep  th 
dep  th-dep th+snpp  tn 
depiiix2  (1,2)  -avs  tp 
depnix2(3 , 2)-avbottp 
depmx2 (3,1) -dep  th 
depmx2(2 ,  l)-thkinx2(l ,  1) 

depmx2  (2 , 2)— (  (depnix2  (2 ,  l)/depth)*(avbottp  -avstp)  )+avscp 
k25-l 

wrlte(9 ,*) ' FLAG  4' 

goto  3000 

else 

if (nomatl.eq.l)then 

if(thkiiixl(l,2)  .eq.vnewsn(l,2))then 
we  1- ( thknocl  ( 1 , 1 )  *10)  *vnewsn  ( 1 , 6 ) 
depChl-(wel/voldsn(l,6))/10 
depth-depthl 

c  write(9,*) 'wel,depthl, depth-  ' .wel . depthl , depth 

end  if 

write(9,*)'FLAG  5' 

goto  200 

else 

we  1- ( thlanx2  ( 1 , 1 )  *10  )  *vnewsn  ( 1 , 6 ) 
depthl-(wel/voldsn(l , 6) )/10 
dep  th- thknix2  (1,1)  +dep  thl 
200  nomatl-1 

nipts-2 
do  250  J-1,6 
thkmxld ,  j  )-vi)ldsn(l ,  j  ) 

250  continue 

thkmxl (1,1) -dep  th 
k25-k25+l 
wrlte(9,*)'FLAG  6' 
goto  4000 
end  if 
end  if 
else 

if  (snpp  tn.gt.  0.0)  Chen 
lf(k2S.eq.l)then 
if (nomatl . eq . 1) then 

wrlte(9,*)'l  thkmxl  (l.D-  ’ ,  thkmxKl,  1) 
write(9,*)'l  snpptn-  ’ ,snpptn 
if(k2.eq.l)then 
thkmxl(l,l)-0.0 
end  if 


c 


c 


300 


350 


chkmxl ( 1 . 1 ) - chkmxl ( 1 , 1 ) +snpp  tn 
write(9.*)'6  depmxl(2,l)-  '  ,  <iepiiixl(2 , 1) 
wrice(9,*)'2  chlanxl(l ,  1)-  '  ,  thkmxKl ,  D 
k25-l 

write(9.*)'FlJ^G  7' 

goto  4000 

else 

thknix2  ( 1 , 1 ) -Chkmic2  ( 1 , 1 ) +snpp  tn 
depth— thknix2  (1 ,  l)+thkiiix2(2 ,1) 
k25-l 

write (9,*) 'FLAG  8’ 
goto  3000 
end  1£ 


else 

if (nomat 1 . eq. 1) then 

we 1- ( thkmxl ( 1 . 1 ) * 10 ) *vnewsn( 1 , 6 ) 
depthl-(wel/voldsn(l , 6) )/10 

depth-dep chi 
wriCe(9 ,*) ' FLAG  9' 
goto  300 
end  if 

vnewsn(l , 6)-0. 280 

we 1- ( thkmx2 ( 1 , 1 ) * 10 ) *vnewsn ( 1 , 6 ) 

vnewsn(l , 6)-0 . 280 

dep  Chi- (we 1/vo Idsn (1,6))/10 

dep  ch-dep  Ch 1+ Chkmx  2(2,1) 

nomaCl“2 


nipCs-3 
do  350  j-1,6 
chkiiix2  ( 1 ,  j  )  -vnewsn  ( 1 ,  j  ) 
chknix2  ( 2 ,  j )  -voldsn(  1 ,  J  ) 


continue 

thkiiK2  (1.1)  -snpp  tn 
thkiiK2  ( 2 , 1 )  -dep  ch+snpp  tn 
depth-thkmx2 (2.1) 
depmx2(l,2)-avstp 

depmx2 (3,2) -avbo  c  cp 
depox2 (3,1) -dep  th 
depnnc2(2,l)-Chkmx2(l.l) 

depnK2(2,2)-((depmx2(2.1)/depth)*(avbotcp-avstp))+avs 

k25-l 


write (9.*) 'FLAG  10' 
goto  3000 
end  if 


else 

if  (nomad .  eq .  1)  then 
k25-k25+l 

write (9.*) 'FLAG  11' 

goto  4000 

else 

k25-k25+l 

write (9 ,*)' flag  12' 

goto  3000 
end  if 


end  if 
end  if 


c 


c**ManipulaCion  of  the  snowpack  to  allow  for  the  effects  of  snowmelt, 
c 

3000  if (depth. le . 5 . 0) then 

write(9,*) 'TERMINATION  DUE  TO  COMPACTION' 

goto  5002 

else 

goto  3001 
end  if 

3001  if(nirsn.ge.depth)then 
write (9,*) 'FLAG  13' 
goto  5002 

else 

if (mrsn. It. 1.0) then 
depnBc2  ( 1 , 2  ) -avs  tp 
depnix2(3 , 2)-avbottp 
depnix2(3 , 1) -depth 
depmx2  (  2 , 1 )  -  thkiiix2  ( 1 . 1 ) 

depnix2(2 , 2)-(  (depmx2(2 ,  l)/depth)*(avbotcp-avstp)  )+avstp 
wriCe(9,*)'N0  MELT,  DEPTH  MATRIX: ' 
do  375  1-1,3 

wrice(9,555)(depiiix2(i,j),j-l,2) 

375  continue 

write (9,*) 'FLAG  14' 

goto  3500 

else 

if  (mrsn.  ge .  thkinx2  ( 1 , 1  > )  then 
nomatl-1 
nipCs-2 
do  400  j-1,6 
thkinxl(l ,  j  )-voldsn(l,  j ) 

400  continue 

dep  th-dep  th - mr sn 

thkmxl (1,1) -dep  th 

depmxl (1,2) -avs  tp 

depmxl(2 , 2)-avbottp 

depmxl  (2 ,  l)-thkiiixl(l ,  1) 

write ( 9 , *) ' ****************  > 

write(9 ,*) ' Post  melt-depth  matrix' 

do  410  1-1,2 

write(9, 555) (depmxl (i,j) ,j-l,2) 

410  continue 

write (9,*) 'FLAG  15' 

goto  4500 

else 

thkmx2 ( 1 , 1 ) -thkmx2 ( 1 , 1 ) -mr sn 

depth-thkmx2(l,l)+thkmx2(2 , 1) 

depmx2 (1,2) -avs  tp 

depaix2  ( 3 , 2  )-avbottp 

depmx2 (3,1) -dep  th 

depmx2 ( 2 , 1) -thkmx2 ( 1 , 1 ) 

depmx2(2, 2)-( (depmx2(2 , l)/depth)*(avbottp-avstp) )+avstp 


415 

3500 


420 


4000 


4001 


425 


430 


wr  1 te ( 9 , *) ' **************** ' 
write(9 ,*)' Post-melt  depth  matrix' 

do  415  1-1,3 

write(9.555)(depmx2(l.J),J-1.2) 
continue 
end  If 

write(9 ,*)' ***************** 
wrlte(9',*) 'THICKNESS  MATRIX' 
do  420  1-1,2 

write(9.550)(thknix2(i.j).j-l,6) 

write(6,550)(thknix2(i,j),j-l,6) 

continue 

write  (9  ,*)  'FLAG  16' 
goto  5000 
end  If 
end  if 

if (depth. le. 5.0) then 

write (9,*) 'TERMINATION  DUE  TO  COMPACTION' 

goto  5002 

else 

goto  4001 
end  if 

i f (mrsn . ge . dep  th ) then 
write (9,*) 'FLAG  17' 
goto  5002 
else 

if (mrsn. It. 1.0) then 
dcpmxl(l,2)-avstp 
depniacl(2 ,2)-avbottp 

depmxl (2,1) -thkmxl (1,1) 
write(9,*)'N0  MELT,  DEPTH  MATRIX:' 
do  425  1-1,2 

wrice(9,555)(depmxl(i,j),j-l,2) 

continue 

write(9 ,*) ' FLAG  18' 

goto  4500 

else 

depth-depth-mrsn 
thkmxl(l , l)-depth 
depmxl (1, 2 )-avstp 
depmxl(2 . 2)-avbottp 
depmxl ( 2 , 1 ) -thkmxl ( 1 . 1 ) 

wiite(9,*) '*******************' 

write(9,*)'Post-melt  matrix;' 

do  430  i-1,2 

write(9, 555)(depntxl(i,  j) ,  j-1,2) 
continue 
end  if 
end  if 


4500 


wrice(9!*VTHICKNESS  MATRIX:' 


write(9.550)thkiiixl(l,l)  .thknocld.Z)  .thlanxl(l,3)  ,thlaBxl(l,4) . 

&  thkmxl(l,5),thkaxl(l,6) 

write(6.550)thkiiixl(l.l).t:hknixl(1.2),thknixl(1.3)  ,thknixl(l,4)  , 

&  thknixl(1.5),thknixld,6) 

5000  write(6,*)'nomatl,nipts-  ' .nomatl .nipts 

write(9,*)'noiiiatl,nlpts-  '  ,nomatl,nipcs 
wrlte(6.*)'k25,kl.k2,k3-  ' .k25.kl,k2,k3 
write(9,*)'k25,kl,k2,k3-  ' .k25 ,kl ,k2 ,k3 
write(6 ,*) 'depth-  '.depth 
write (9.*) 'depth-  '.depth 

wrlte(6,*) 'snowfall. meltrate(cni  snow)-  '  . snpptn.mrsn 
wrlte(9.*)'snowfall,meltrate(cm  snow)-  '.snpptn.mrsn 
write(9.*) 'avstp.avbottp-  ’ .avstp . avbottp 
550  foniiat(6f6 . 2) 

555  format(2f6 . 2) 

c 

c**Read  in  next  dally  data  line,  shift  existing  bottom  data  line  to 
c  the  top . 
c 

9000  yyy(2.6)-yyy(2.6)/100.0 

yyy (2 . 2 )-yyy ( 2 . 2 ) *100 

cuintin-xxx(2 , 1) 

write(9,*) '******cumulative  time-  '.cumtm 

xxx(l,8)-xxx(2.8) 

xxx(l . l)-time2 

yyy<i-i)"yyy(2,i) 

yyy(1.2)-yyy(2,2) 

yyy.'1.9)-yyy(2.9) 

yyy(l.3)-yyy<2,3) 

)07(i'6)"yyy(2.6) 

yyy(1.7)-300^(2,7) 

xxx(2.8)-0 

xxx(2 . l)-0 

yyy(2,l)-o 

yyy(2.9)-0 

yyy(2.2)-0 

yyy(2,3)-0 

yyy(2.6)-0 

yyy(2,7)-0 

nkont-n-kontrol 

if (nkont . eq . 1 ) then 

goto  5002 

end  if 

read(10.*)xxx(2,8),xxx(2.1).yyy(2.1) ,yyy(2,9) .yyy(2.2) , 
yyy(2,3) , 
yyy(2.6).yyy(2.7) 


& 

& 


88 


5001 

5002 


wrlte(9 ,*)' shifted  matrix  to  read  next  line  ’ 
write(6,88)xxx(l,8) ,xxx(l,l) ,yyy(l,l) ,yyy(l,9) ,yyy(l,2) , 
&  yyy(l,3), 

yyy(i.6)  .yyy(1.7) 

write(9,88)xxx(l,8) ,xxx(l,l) ,yyy(l,l) ,yyy(l,9) ,yyy(l,2) , 
&  yyy(i.3), 

yyy(1.6),yyy(l,7) 

write(6,88)xxx(2,8) ,xxx(2,l) ,y3o^(2, 1) ,yyy{2 ,9) ,yyy (2,2) , 
&  yyy(2.3), 

S'  yyy(2.6>  .yyy(2.7) 

write(9 ,88)xxx(2,8) ,xxx(2,l) ,yyy(2,l) ,yyy(2 ,9) ,yyy(2 , 2) , 
S'  yyy(2,3). 

S'  yyy(2,6)  ,yyy(2,7) 

format(7f6.1) 
airtp-yyy(l,l) 
yyy(I.l)-y-yy(-.l)-273.l5 
yyy(i.9)-yyy(3.9)-273.i5 
kelvln-yyy(l,l)+273.15 
PPtn-yyy(l,7) 
ddate-xxx(l,8) 
kontrol-kontrol+1 
snpptn-0 .0 

write(9,*)'CALL  NEXT  L' 
write(9,*)'  * 
write(9 ,*) 
write(9 .*) 
continue 
write(9 ,*) 
stop 
end 


FLAG  30' 


PROGRAM  TERMINATION' 


********** 


C 


c***'*‘*'**subroutine  cmelt******** 

c  This  calculates  Che  melcrace  In  cm  of  snow  and  mm  of  water 
equivalent. 

c  The  density  of  the  snowpack  is  calculated  if  the  pack  is  a  2- 
layered  pack. 

c  Melt  is  calculated  using  Che  snow  energy  budget.  The  components  of 
the 

c  snow  energy  budget  are  calculated  by  TSTM,  rainmp  and  dU/dtis 
calculated 
c  within  cmelt. 
c 

subroutine  cmel t(l,n,enavm,ne trad, p term, tsol , tgbr , tabsor , 

(t  taterm,  chterm,  cdterm,pptn,densn,mrwe  ,mrsn,  gterm,  denw,  depth , 

&  avsCp,avboCCp) 

conimon/matrlx2/depnixl(2,2)  ,depmx2(3 , 2) 

common/maCrlx3/Chkmxl(l ,6) , thkmx2(2 , 6) 

common/tsCm2/nomatl 

common/phy c  ons /lheatf,hcapsn,hcapw 

real  Iheatf ,ne trad, mrwe ,mrsn, depth, densn, A, B, C , D, avsntp , dudt 
c 

c**calculaClon  of  the  density  of  the  snow, 
if  (nomad .  eq .  1)  Chen 
densn-thkmxl (1,6) 
else 

A-thkmx2(l, 1) /depth 
B-Chkmx2(2 , l)/depth 
C-A*chkmx2(l,6) 

D-B*thkDix2(2,6) 

c  wrice(9,*) 'A,B,C,D,depth-  ' ,A,B,C,D,depCh 

c  wrlte(9,*) 'thkmx2(l,l)  ,thkmx2(2,l)“  '  ,  thknix2(  1 , 1) ,  Chkmx2 (2 , 1) 

c  wrlte(9,*) ' chkmx2(l ,6) , thkmx2(2 , 6)-  ' , thkmx2(l , 6) , thkmx2 (2 , 6) 

densn-C+D 
end  if 

wrlte(9,*)'  ' 

write(9,*)'l  density  of  snow-  ', densn 
c 

c**Calculation  of  dU/dt,  J/m**2.day 
avsntp- (avstp+avbottp)/2 

dudc-depth*( ( (densn*1000)*hcapsn)+hcapw)*avsntp 
write(9 ,*) ' dU/dt-  ’ ,dudt 
c 

c**Calculation  of  melcrace. 

netrad-tsol- (tsol-tabsor)+taCerm- tgbr 
enavm-netrad+thterm+tdterm+gterm 
enavm-enavm*8 . 64e+4 
c  lW/m**2  -  8 . 64e+4  J /m**2 . day 

wrlte(9,*) 'ENAVM  1-  ' ,enavm 

wrlte(9 ,*) 'TSOL,TABSOR,TATERM,TGBR-  ' , tsol , tabsor , taterm, tgbr 
write(9 ,*) 'NETRAD,THTERM,TDTERM-  ' ,ne trad, chterm, tdcerm 
write(9,*) 'GTERM, PTERM-  ’, gterm, pterm 
enava-enavni+pcerm-dudc 
write (9,*) 'ENAVM  2-  ' .enavm 
mrwe- ( enavm/ ( Ihea tf*denw) ) *1000 
mr  sn-  (  enavm/  ( IheaC  f  *  (  densTi*1000 ) ) )  *100 
c  Multiply  mrwe  by  1000,  so  that  result  is  in  mm  water, 
c  Multiply  mrsn  by  100,  so  chat  result  Is  In  cm  snow, 
return 
end 


c****‘*'**subroutine  clmelt******** 

c  This  is  identical  to  cmelC  except  that  daily  averages  are  used  in 
the 

c  calculation  of  melt  not  totals, 
c 

subroutine 

clmeltd  ,n,  enavml ,  ne trad, p term, avsol, avgbr ,  avabsor , 

& 

avaterm , avhterm , avdterm , pptn , densn, mrlwe , mrlsn , gcerm , denw , depth , 

&  avstp , avbottp) 

coiiimon/maCrix2/depmxl(2 ,2)  ,depmx2(3 ,2) 
common/matrix3/thkmxl(l,6) , thkmx2(2, 6) 
coomon/tstm2 /nomad 
common/phycons/lheatf ,hcapsn,hcapw 

real  Iheatf ,ne trad, mrlvfe , mrlsn, depth, A, B, C, D , avsntp, dudt 
c 

if  (nomad  .  eq .  1)  then 
densn-thkmxl (1,6) 
else 

A-chkmx2(l, 1) /depth 
B-Chkmx2(2 , l)/depth 
C-A*thkmx2 (1,6) 

D-B*thkmx2(2,6) 

write(9 ,*) ' A, B,C , D, depth-  ' .A, B,C , D , depth  " 

write(9,*) '  thlonx2(l,  1)  ,  chlaiix2(2 , 1)-  '  ,  thlanx2(l ,  1)  ,  thlanx2(2 . 1) 

write (9,*) ' chkmx2(l,6) , thkmx2(2 , 6)-  ' , chkmx2(l , 6) , thlanx2(2 ,6) 

densn-C+D 

end  if 

write(9,*)'  ' 

write(9,*)'2  densn-  ', densn 
c 

avsntp- ( avs  tp+avbo  t tp ) /2 

dudt-depth*( ( (densn* 1000)*hcapsn)+hcapw)*avsntp 
wrice(9,*) 'dU/dt-  ',dudc 
c 

ne trad-avs o 1 - ( avs o 1 - avab s or ) +ava te rm - avgb r 
enavml-necrad+avatarm+avdterm+gterm 
enavml-enavml*8 . 64e+4 
wrlte(9 ,*) ' ENAVM  3-  ', enavml 

wrice(9,*) 'AVSOL, AVABSOR, AVATERM-  ' , avsol , avabsor , avaterm 
wr ite (9 ,*) 'NETRAD, AVHTERM, AVDTERM-  ' ,ne trad, avhterm, avdterm 
writa(9,*)'GTERM,PTERM, AVGBR-  ' ,gterm,pterm, avgbr 
enavml-enavml*p tera-dudc 
wrlte(9,*)' ENAVM  4-  ', enavml 
mrlwe-(enavml/(lheatf*denw) )*1000 
mrlsn-(enavml/(lheatf*(densn*1000)*0. 96) )*100 
return 
end 
c 
c 


C*******************************************************************^ 

*★*★****★ 

c 

c***-****subroucine  ralnaip******* 

c  This  calculates  the  energy  input  to  the  pack  by  rain  falling  on 
c  a  melting  pack  at  0  oc. 
c 

subroutine  rainmp(i , n.pptn.pterm, denu) 
common/rain/airtp 

common/phycons/lheatf , hcapsn , hcapw 
real  pterm 

pCenB-(denw*hcapw*aircp*ppcn)/1000 

return 

end 

c 

c 


******* 


c*******subroutine  rainfp******* 

c  This  calculates  the  warming  effect  of  rainfall  on  a  snowpack  with 
c  a  temperature  <0  oC. 

c  Iheatf  converted  from  HJ/kg  to  kJ/kg,  rstpc  calculates  in  kJ/kg 
c  densn  converted  from  g/m**2  to  kg/m**2 
c  depth  converted  from  cm  to  m 
c 

subroutine  ralnfp(i,n,pptn, rstpc , denw, avstp , depth , densn) 
common/rain/air tp 

common/phycons/lheatf , hcapsn . hcapw 
common/matrix2/depmxl(2 ,2) ,depmx2(3 , 2) 
common/matr ix3/thkmxl ( 1,6), thkmx2 (2,6) 
common/tstm2/nomatl 

real  rstpc, rstpcl , rstpc2 , rstpc3 , Iheatf , depth, densn 

depth-depth/100 

densn-densn*1000 

lheatf-333 . k 

rstpcl-(hcapw*airtp)+lheatf- (hcapsn*avstp) 

rstpc2-hcapsn*densn*depth 

rstpc3-( ( (pptn/1000)*denw)*rstpcl)/rstpc2 

rstpc-rstpc3 

avs  tp-avs  tp+r s  tpc 

c  write(9 ,*)' rstpcl , rstpc2 , rstpc3 , avstp-  ', rstpcl , rscpc2 , 

c  &  rstpc3, avstp 

lheatf-0 . 3  34e+6 
densn-densn/1000 
depth-depth*100 
return 
end 
c 
c 

★★★*★★★★★★*★★**★★**★★*★★★★★**★★*★★***★★★★★★★★★■**’**★★■**★•**■**'*★★*■*'***** 

******* 


c 


c*****'**subroutine  csfrq******* 

c  This  calculates  sfrq.  Sfrq(ix)  is  the  vertical  grid  spacing  in 
c  cm  in  each  layer  tx  in  cm**2/m.  A  minimum  value  for  sfrq  of 
c  5cm  is  used.  Numerical  instability  results  if  smaller  values 
c  are  used, 
c 

subroutine  csfrq(i,n) 
common/matrix3/thkmxl(l ,  6) ,  thkji!x2(2 , 6) 
c ommon/ t s tm2 /noma t 1 

common/vegi/iveg, sigf , state , epf , fola ,hfol 
common/s  frq/s  f rq 

common/ ts tm3/to tt im,  tfrq ,  tpmt  .nipts 
real  sfrq(6) , sfrql . sfrq2 
if  (nomad .  eq .  1)  then 
sfrql-(Chkiiixl(l.l))/2 
else 

sfrql-(thknix2(l,l))/2 
end  if 

100  if(sfrql.le.5.0)then 

goto  200 
else 

if (sfrql . le . 10 . 0) then 

goto  200 

else 

sfrql-sfrql/2 
goto  100 
end  if 
end  if 

200  sfrq(l)-sfrql 

if(nomaCl.eq. 2) then 
sfrq2“thknix2(2 ,  l)/2 
300  if(sfrq2.1e.5.0)then 

goto  400 
else 

if(sfrq2.1e.l0. 0) then 

goto  400 

else 

sfrq2-sfrq2/2 
goto  300 
end  if 
end  if 

400  sfrq(2)-sfrq2 
end  if 

timespl-0. 5*(sfrq(l)*sfrq(l)/chkraxl(l ,2) ) 
print  *,'tinie  '  ,  timespl ,  sfrq(l) ,  thkmxl  (1 . 2) 
if (nomatl . gt . 1)  Chen 

timespl-0. 5*  (sfrq  (l)*sfrq(  l)/thlaiix2(l ,  2) ) 
tlmesp2-0. 5*(sfrq(2)*sfrq(2)/thlcmx2(2 . 2) ) 
timespl-min( timespl , tlmesp2) 


end  if 


itime-lnt( timespl) 
print  ♦,'lcime  '.itlme 
if ( itime . eq . 0) itime-1 
if(itime.gt.  5)itinie-5 
if (iveg. eq. 1) itime-1 
tfrq-real( itime) 
return 
end 
c 

'k'k'k'k'k'k'k’k'k 


c******Subroucine  TSTM******* 

c  This  calculates  the  longwave,  shortwave,  latent  heat  and  sensible 
c  heat  exchanges  over  the  snowpacVc  or  soil  surface,  and  the  heat 
c  transport  through  the  snowpack  or  soil  pack.  TSTM  has  been 
converted 

c  to  a  subroutine  from  BALICk,  L.K. ,LINK,  L.E.,  SCOGGINS,  R.K.  and 
c  SOLOMAN, J , 1 . (1981)Thermal  Modelling  of  Terrain  Surface 
Elements . ,U. S . 

c  Army,  Tech.  Report  No.  EL-81-2,  Washington  D.C.  and  BALICK,  L.K. , 
c  SCOGGINS,  R.K.  and  LINK,  L.E. (1981)Inclusion  of  a  simple  vegetation 
c  layer  in  terrain  temperature  models  for  thermal  infrared(IR) 
c  signature  prediction. , Miscellaneous  Paper  EL-81-4,  U.S.  Army 
c  Engineer  Waterways  Experiment  Station, CE,Vicksberg, Miss .  Conversion 
c  to  a  subroutine  was  facilitated  by  the  help  of  Mr.  Randy 
Scoggins (WES) . 
c 

SUBROUTINE  TSTM( 1 ,anl , an,klO,avstp , avgbr , avsol , avabsor , 

&  avaterm, avhterm, avdterm,k2 .avbottp , tgbr , tsol , 

&  tabsor , taterm, thterm, tdterm) 

C  TSTM:VEGIE .  BARE  OR  VEGETATED  SURFACE  TEMPERATURE  MODEL 

C . 

c- 

C  THE  FUNCTION  OF  THIS  PROGRAM  IS  TO  PREDICT  SURFACE 

C  TEMPERATURE . 

C 

C . 

C  PARAMETER  (ND-30) 

INTEGER  OUTGO 
DIMENSION  datmx(100,10) 

DIMENSION  RHOC(6),MAX(10),DEPTH(450) ,PROF(2,450) , 

&  FMM(30,10),BBB(30,10) 

DIMENSION  TITLE(7) 

DIMENSION  CLR(8) ,NX(6) ,ATF(2) ,FEB(2) 

DIMENSION  RR(6) ,INTR(7) 

dimension  chk(6) , alph(6) sCor(7 ,450) , stabn( 6) 

REAL  KTEMPG,KTEMPA,LAT,ACL(8) ,BCL(8) ,M,KSQ,bottp 
CHARACTER  DATE*8 ,HEA0ER*72 , AN*1 , AN1*1 
CHARACTER*30  FNAME 


conimon/matrixl/xxx(30 , 10)  ,yyy(30 , 10) 
common/matr  lx2/depmxl  (2,2),  depaix2  (3,2) 
conunon/matr  ix3/thkiiixl  (1,6),  chla]ix2  (2,6) 
common/tstml/press , ncloud, za , slope , surfac , day , laC , vec 
common/ts tm2/nomatl 
common/tstmS/tottlm,  tfrq ,  tpmt  ,nipts 
connnon/vegl/iveg ,  s  Igf ,  s  taCe ,  cpf ,  f  ola ,  hf  ol 
conmon/sfrq/sf rq (6 ) 

DATA  CLR/0. 04,0.08,0.17, 0.20, 0.22, 0.24, 0.24, 0.25/ 

DATA  OUTCD/0/ 

DATA  ACL/82.2,87.1,52.5,39.0,34.7,23.8,11.2,15.4/ 

DATA  BCL/.079, . 148, . 112, .063, .104, .159, - .167, .028/ 
DATA  SIGMA, PI , AC,BC/8 . 12E-11, 3 . 141593 , 17 . 269 , 35 . 86/ 
DATA  CC/0.261/ 

DATA  LAST, G,KSQ, CP/2, 980. 0,0. 16, 0.24/ 

C  SET-SPECIFICATIONS 

C  STATEMENT  FUNCTIONS  FOR  USE  IN  VEGETATION  SECTION 

E(T)-RH*6 . 108*EXP(AC*(T-273 . 15)/(T-BC) ) 

ESAT(T)-6 . 108*EXP(AC*(T-273 . 15)/(T-BC) ) 

Q(T)-0. 622/(PRESS/E(T) - . 378) 

QSAT(T)-0 . 622/(PRESS/ESAT(T) - . 378) 
print  Inside  tstm' 

SOLCHK-0 
lflag2-0 
lrecm-0 
ipmt-0 
sun-0 . 0 

write (9,*) 'data  matrix  received  by  tstm  ' 
write(9 ,60)xxx(l, 8) ,xxx(l,l) ,yyy(l,l) ,yyy(l , 9) ,yyy(l , 2) 
&  yyy(i.3). 

&  yyy(l,6) ,yyy(l,7) 

write(9,60)xxx(2,8) ,xxx(2,l) ,yyy(2,l) ,yyy(1.9) ,yyy(2 ,2) 
&  yyy(2,3), 

yyy(2,6)  ,yyy(2,7) 

60  format(7f6 . 1) 

if (nomatl . eq. 1) then 
xxx(l ,  5)-depnixl(l ,  1) 
yyy(l,5)-depmxl(l,2) 
xxx(2 , 5)-depmxl(2 , 1) 
yyy(2,5)-depmxl(2 ,2) 
thk( 1 ) -thkmxl ( 1 , 1 ) 
alph(l)-thknixl(1.2) 
fk(l)-thkmxl(l, 3) 
c  write(9,*)'sfrq  1  ,sfrq 

epsn- thkmxl ( 1 , 4 ) 
smalla-(l- thkmxl (1 , 5)) 
c  write(9 ,*) 'TPRNT-  ’ , tprnt 

c  wrlte(9,*)'ABS0RPTIVITY-  ' ,smalla 

c*^smalla-l-albedo(see  Balick  et.  al.  original) 
else 


XXX  ( 1 , 5 ) -depinx2  ( 1 , 1 ) 
yyy(l ,  5)-depnix2(l ,  2) 
xxx(2,5)-depiiix2(2,l) 
yyy(2,5)-depinx2(2,2) 
xxx(3 , 5  )-depiitx2  ( 3 , 1 ) 
yyy(3,5)-depnix2(3,2) 

thk(l)-thkiiix2(l.l) 

thk(2)-t:hkmx2(2.1) 

alph(l)-thkiiix2(l,2) 
alph(2)-thkiiix2(2,2) 
fk(l)-thkmx2(l,3) 
fk(2)-thkiiix2(2,3) 
write(9 ,*) ' sfrq-  2',sfrq 
epsn-thknix2  ( 1 ,  A ) 
smalla- ( 1 - thkmx2 (1,5)) 

write (9,*) '**2  lAYER  SYSTEM  INPUT  INTO  TSTM**' 
c  write(9,*) 'xxx(l,5) ,yyy(l,5)-  ' ,xxx(l , 5) ,yyy (1 , 5) 

c  write(9,*)'xxx(2,5) ,yyy(2,5)-  ' ,xxx(2 , 5) ,yyy(2 , 5) 

c  write(9,*)'xxx(3,5) ,yyy(3,5)-  ' .xxx(3 , 5) ,yyy(3 , 5) 

c  wrlte(9 ,*) ' thk(l) , chk(2) , alph(l)-  ' , thk(l) , thk(2) , alph(l) 

c  write(9 ,*) ' alph(2) ,fk(l),fk(2) ,epsn- 

' .alph(2).fk(l),fk(2),epsn 
c  write(9,*) 'TPRNT-  '  ,  tpmt 

c  write(9,*)'ABSORBTIVITY-  ' .smalla 

end  if 

C . 

C  INITIALIZE-VARIABLES -AND- CONSTANTS 

99996  ASSIGN  99994  TO  199995 

GO  TO  99995 
99994  CONTINUE 

C . 

C  INPUT -DATA 

ASSIGN  99989  TO  199990 
GO  TO  99990 

C . 

C  PRINT -INPUT -DATA 

99989  ASSIGN  99987  TO  199988 

IF(AN1.EQ. 'Y')GO  TO  99988 
99987  CONTINUE 

C . 

C  CALCULATE-TABLE-SLOPE-AND- INTERCEPT 

99986  ASSIGN  99982  TO  199983 

GO  TO  99983 

C . 

C  CALCULATE-SURFACE-AND-LAYER-TEMPERATURE 

99982  CONTINUE 

ASSIGN  99980  TO  199981 
GO  TO  99981 

C . 


no  on 


99980  CONTINUE 

PPIN-C-mCI ,  7)+YYY(2 , 7)  )/2 

sunic2-0 

sumc3“0 

sunc4-0 

suiiic5“0 

sumcS-O 

sumc7-0 

suncS-O 

suiiic9-0 

do  99  1-1, klO 

sumc2-sxiinc2-t-datmx(l ,  2) 

sunic3-sunic3+datnix(  i ,  3 ) 

sumc4-suiiic4+<Ucmx  ( 1,4) 

suiiic5-sumc5+dacmx(  1 , 5 ) 

sumc6-suiiic6+datmx  ( 1 , 6 ) 

sumc7-sumc7+daCmx( 1 ,7) 

sumc8-suiiic8+datmx(  1,8) 

suinc9-suinc9+datnix(  1 ,9) 

99  conclnue 

Cgbr-sumc3 

Csol-sunic4 

tabsor-sumcS 

tatem-simcS 

Chtenn-sunic7 

Cdterm-sunic8 

avscp"sunic2/lcl0 

avgbr-suiiic3/kl0 

avsol-suiiic4/kl0 

avabsor-sumc5/kl0 

avacerm-suiiic6/kl0 

avhtejrm-sumc7/kl0 

avdtertn-sumc8/kl0 

avboctp-suiiic9/kl0 

WRITE(*,*) 'NORMAL  TERMINATION' 


STOP 

RETURN 


99997  CONTINUE 

C  FORMATS 

90  FORMATC'  BOTTOM  BOUNDRY  INDEX-', 13) 

92  FORMATC'  BOTTOM  BOUNDRY  TEMPERATURE- ', F6 . 1 , '  DEG  C’) 
95  FORMATC'  BOTTOM  BOUNDRY  HEAT  FLUX-' , F6 . 1 , 'W/M**2 ' ) 

97  FORMATCFll . 2 , Fll . 1 , F13 . 2 , F13 . 1 , F8 . 2) 

120  FORMATCA8,F6,l) 

140  FORMATCF5 , 1 , FIO . 1 , F6 . 1 , F12 . 1 , F12 . 1 . F13 . 1) 

145  FORMATCF6.1,F7.1) 

150  F0RMATCF12.1,I12,F11.2) 

160  FORMATCFll, 2, F13. 2, F8.1, FIO. 1) 

170  FORMATCF9.4,F12.1) 

180  FORMATCI7,F16.1,F11.0,F12.0) 

190  F0RMATC4X,F5.2,F6.1,2X,I8,I9,I11.I9,I9,I8) 

260  FORMAT  C  9X , F8 . 1 , FI 1 . 2 , F9 . 2 , F9 . 2 , F8 . 2 , FIO . 2 , F8 . 2 ) 

200  F0RMATCF6.2,F8.2,F17.2) 

210  F0RMATCI6,F11.1,F12.1,F14.2,F14.2) 

220  FORMAT ClHl) 


230  FORMAT(A72) 

240  FORMAT(4X,F7.2,3X,F5.0,5X.F6.2.10X,F6.2,9X,F7.2) 

250  FORMAT(6X,F6.1,13X,F4.1.12X,F5.1) 

310  FORMAT (IIX, 'TOTAL  GRAYBODY  EFFECTIVE  GROUND  FOLIAGE', 

&  4X, 'SOLAR') 

320  FORMAT (14X, 'RADIANCE  TEMP' . lOX, 'TEMP  TEMP' 

&  .4X, 'INSOLATION') 

330  F0RMAT(5X. 'HR' .7X. ' (W/M**2) ' ,9X, ' (C) ' ,11X, ' (C) ' 

&  .6X,'(C)  (W/M**2)') 

340  F0RMAT(9X.  ' . REFL-NREFL - REFL - NREFL'  ,30(1H-)) 

270 

F0RMAT(3X. F5 . 2 , 5X, 14 , 2X, 14 , 3X, F6 . 1 . 2X, F6 . 1 , 3X, F6 . 1 , 3X, F6 . 1 , 7X,  14) 

350  FORMAT(57X, 'SENSIBLE  LATENT') 

360  F0RMAT(5X, 'HRS  SURFACE  GRAYBODY  SOLAR  SURFACE  ATMOS', 
&  '  IR  HEAT  HEAT') 

370  FORMAT (1 IX, 'TEMP  RADIANCE  INSOLATION  ABSORP  EMISSION', 

&  '  LOSS  LOSS') 

380  F0RMAT(11X. 'DEG  C' ,23(1H- ) , ' (M/M**2) ' , 24(1H- )  ) 

390  FORMAT (3X, 'TIME' ,18X, 'DEPTH,  TEMPERATURE ' /4X ,' HR ',  2 IX ,' CM '  , 

&  9X, 'C'/2X,65(1H-)) 

400  F0RMAT(2H0  , F5 . 2 ,4(3X,F5 . 1 , '  ',F5.2)) 

410  FORMAT( lOX , F5 . 1 , IX , F5 . 2 , 3X , F5 . 1 , IX, F5 . 2 , 3X , F5 . 1 , IX , F5 . 2 , 

&  3X,F5.1,1X,F5.2) 

GO  TO  199997 

C . 

99995  CONTINUE 

C  TO  INITIALIZE-VARIABLES-AND-CONSTANTS 
BB— 2.4E-4 
MAX(l)-2 
MAX(2)-2 
MAX(3)-2 
MAX(4)-2 
MAX(5)-NIPTS 
MAX(6)-2 
MAX(7)-2 
MAX(8)-0 
MAX(9)-0 
MAX(10)-3 
IBUG-0 
IEOF-0 

GO  TO  199995 


99990  CONTINUE 
C  TO  INPUT -DATA 
C 

C  INPUT-ATMOSPHERIC-SPECIFICATIONS 

ASSIGN  99978  TO  199979 
GO  TO  99979 
C 

C  INPUT-SURFACE-ORIENTATION-SPECIFICATIONS 

99978  ASSIGN  99976  TO  199977 
GO  TO  99977 

C 

C  INPUT-HEAT-FLOW-CACULATION-CONTROLS 

99976  ASSIGN  99974  TO  199975 
GO  TO  99975 
C 

C  INPUT-INITIAL-TEMPERATURE-PROFILE 

99974  ASSIGN  99972  TO  199973 
GO  TO  99973 
C 

C  INPUT -TOP- SURFACE -CONSTSNTS 

99972  ASSIGN  99970  TO  199971 
GO  TO  99971 
C 

C  INPUT-LAYER-SPECIFICATIONS 

99970  ASSIGN  99968  TO  199969 
GO  TO  99969 
C 

C  INPUT-BOTTOM-BOUNDRY-DATA 

99968  ASSIGN  99966  TO  199967 
GO  TO  99967 
C 

C  INPUT- VEGETATION- PARAMETERS 

99966  ASSIGN  99798  TO  199799 
GO  TO  99799 
C 

99798  GO  TO  199990 

C . 

99979  CONTINUE 

C  TO  INPUT-ATMOSPHERIC- SPECIFICATIONS 

C 

C  LINE  TIME  AIR  TEMP  RH  CLOUD  COVER  WIND  SPEED , INSOLATION 

C  NO.  tai  DEG  C  %  (0-1)  M/3  CAL/CM**2- 

MIN 
C 

xxx(2, l)-xxx(2 , l)+2400*(xxx(2. 8) -xxx(l , 8) ) 

860  DO  99965  J-1,2 

ixxx-lnc(xxx(j ,1)/100.0)*100 
axxx-(xxx(j , 1)  -  real(ixxx))/60. 0 
xxx(j ,l)-real(ixxx)/100.0+axxx 
XXX(J,7)-XXX(J,1) 

XXX(J.2)-XXX(J,1) 

XXX(J.3)-XXX(J,1) 

XXX(J.4)-XXX(J,1) 

XXX(J,6)-XXX(J,1) 

YYy(J , 1)-VYY(J . l)+273 . 1 
YYY(J , 9)-YYY(J , 9)+273 . 15 
YYY(J , 2)-YYY(J , 2)*0 . 01 
YYY(J,6)-YYY(J.6)*100.0 


YYY(J,4)-YYY(J,4)/697.6 
99965  CONTINUE 

XXX(3,10)-XXX(2.1) 

XXX(2,10)-14.0 

XXX(1,10)-XXX(1.1) 

•m(3,10)-YYY(2,l) 

YYY(2,10)-YYY(1,9) 

YYY(1,10)-YYY(1,1) 

PRINT*,  '10s  '  .XXXd.lO)  ,XXX(2,10)  ,XXX(3,10)  .YYY(1,10)  , 

&  YYY(2.10) .YYY(3,10) 

C 

840  FACTH-(1000.0/PRESS)**0.286 

GO  TO  199979 

C . 

99977  CONTINUE 

C  TO  INPUT-SURFACE-ORIENTATION-SPECIFICATIONS 

C 

C  LINE  SFC  SLOPE  SFC  AZIMUTH  DAY  LATITUDE 

C  NO.  DEG-HORIZ-0  DEG  S-O  JULIAN  DEG 

C 

SL0PE-SL0PE*PI/180 . 0 
SURFAC-SURFAC*PI/180 . 

GO  TO  199977 

C . : . 

99975  CONTINUE 

C  TO  INPUT-HEAT-FLOW-CACULATION-CONTROLS 

C 

C  LINE  NO.  OF  NO.  OF  24  HRS  TIME  STEP  P’^TNT  FRE  Q 

C  NO.  LAYERS  REPITIONS  MIN  MIN 

C  1-6  2-5 

C 

TOTTIM-XXX(2,l)  -  XXX(l.l) 

GO  TO  199975 

C . 

99973  CONTINUE 

C  TO  INPirr- INITIAL-TEMPERATURE- PROFILE 
C 

C  LINE  NO.  OF 

C  NO.  PROFILE  POINTS 

C 
C 

C  LINE  DEPTH  TEMP 

C  NO.  CM  DEG  C 

C 

MAX(5)-NIPTS 
DO  99964  J3-1.NIPTS 
YYY(J3,5)-YYY(J3,5)+273.15 
99964  CONTINUE 

GO  TO  199973 

C . 

99971  CONTINUE 

C  TO  INPUT -TOP -SURFACE -CONSTANTS 
C  LINE  EMISSIVITY  ABSORBTIVITY  MOISTURE 

C  NO.  CONTENT  DRY  WT. 

C 

FACTA-SIGMA*EPSN 

C 

GO  TO  199971 


C . 

99967  CONTINUE 

C  TO  INPUT-BOTTOM-BOUNDRY-DATA 
C 

C  IF  LFLUXY-0,  THERE  IS  NO  HEAT  FLUX  THRU  BOTTOM 

C  IF  IFLUXY  LT  0.  THERE  IS  NO  AIRSPACE  BENEATH  BOTTOM 

C  IF  IFLUXY  GT  0,  THERE  IS  AIRSPACE  BENEATH  BOTTOM 

C 

LFLUXY— 1 

IF(.NOT. (LFLUXY.EQ.O))  GO  TO  99962 
DPRM0-YYY(NIPTS.5) 

TB-0. 

GO  TO  99963 

99962  IF(. NOT. (LFLUXY. LT.O))  GO  TO  99961 
C 

DPRM1--0.5 

TB-FK(NOMATL) 

BEP-O . 0 
BK-0. 

REP-0. 

TR-0.0 
FACTD-O . 

FACTE-0 . 

RK-0. 

DPRMl-DPRMl/697,6 
GO  TO  99963 
99961  CONTINUE 

99963  GO  TO  199967 

C . 

99969  CONTINUE 

C  TO  INPUT-LAYER-SPECIFICATIONS 
C 

C  LINE  THICKNESS  VERT.  GRID  THERMAL  DIFF  HEAT  COND 

C  NO,  CM  SPACE-CM  CM**2/MIN  CAL/MIN-CM-K 

C 

DO  99960  J4-1,N0MATL 
RH0C(J4)-FK(J4)/ALPH(J4) 

99960  CONTINUE 

- CONTROL  SWITHCHES  SVECIFIED" 

NIT-1 
IEFSWT-1 
GO  TO  199969 

C . 


99799  CONTINUE 

C  TO  INPUT-VEGITATION-PARAMETERS 
if(iveg.eq.O)  go  to  1120 

print  ' iveg.sigf ,epf .fola.hfol  ' , iveg, slgf , epf , fola ,hfol 
IF(SIGF.LE.O.O)GO  TO  199799 
TF-YYY(1,1) 

IVEG-1 

EPl-EPF+EPSN- EPF*EPSN 
ZO-O . 131*HFOL**0 .997 
CHO-KSQ/ ( ALOG ( ZA/ZO )**2 ) 

ZDSP-0 . 701*HFOL**0 . 979 
CHH-KSQ/(ALOG( (ZA-ZDSP)/Z0)**2) 

CHG-(1. -SIGF)*CHO+SIGF*CHH 
DELTMP-1 . 

QAF-QSAT(TF) 

1120  GO  TO  199799 

C . 

99988  CONTINUE 
C  TO  PRINT-INPUT-DATA 
c  WRITE(*,139) 

WRITE(*,220) 

C  WRITE(*,120)DATE.TTIME 

WRITE(*.*)'  ' 

WRITE(*,230)HEADER 
WRITE(*.*)'  ' 

WRITE(*,*) '  ATMOSPHERIC-SPECIFICATIONS' 

WRITE(*,*)'  ‘ 

WRITE(*,*)'  ATMOS  PRESS  CLOUD  TYPE  SHELTER' 

WRITE(*,*)'  MB  INDEX  HEIGHT-CM' 

WRITE (*, 1 50 ) PRES S , NCLOUD . ZA 
WRITE(*,*)'  ' 

WRITE(*,*)'  TIME  AIR  TEMP  RH  CLOUD  COVER  WIND  SPEED' 

&  ,  '  SOLAR  IRRAD' 

WRITE(*,*)'  HRS  DEG  C  %  (0-1)  M/S 

W/M**2 ' 

MMAX-2 

WRITE(*, 140) (XXX(J . 1) .YYY(J . 1) -273 . 15 . 

& 

YYY(J,2)*100.0.YYY(J,3) , YYY(J . 6)*0 .01, YYY(J , 4)*697 . 6 , J-1 ,MMAX) 
WRITE(*.*)'  ' 

WRITE(*,*)'  SURFACE- ORIENTATION -SPECIFICATIONS' 
WRITE(*,*)'  ' 

WRITE(*,*)'  SFC  SLOPE  SFC  AZIMUTH  DAY  LATITUDE' 

WRITE(*,*)'  DEG-HORIZ-0  DEG  S-0  JULIAN  DEG' 

WRITE(* . 160)SLOPE*180/PI , SURFAC*180 . 0/PI , DAY , LAT 
WRITE(*.220) 

WRITE(*,*) '  HEAT-FLOW-CACULATION-CONTROLS' 

WRITE(*,*)'  ' 

WRITE(*.*)'  NO.  OF  NO.  OF  24  HRS  TIME  STEP  PRINT  FREQ' 
WRITE(*,*)'  LAYERS  REPETITIONS  MIN  MIN' 

WRITE(*,*)'  1-6' 

WRITE ( * , 1 80 ) NOMATL , TOTTIM/24 . 0 . TFRQ , TPRNT 
WRITE(*.*)'  ' 

WRITE(*,*) '  INITIAL-TEMPERATURE- PROFILE' 

WRITE(*,*)'  ' 

WRITE(*, 145) (XXX(J , 5) ,YYY(J . 5) -273 . 15 , 

&  J-1,NIPTS) 

WRITE(*,*)'  ' 


WRITE(*,*)’  TOP-SURFACE-CONSTANTS’ 

WRITE(*.*)'  ' 

WRITE(*.*)'  EMISS  ABSORB  SATURATION' 

WRITE(*.200)EPSN, SMALLA , WET 
WRITEC*,*)'  ' 

WRITE(* . *) '  INPUT-LAYER-SPECIFICATIONS ' 

WRITEC*,*)’  LAYER  THICKNESS  VERT.  GRID  THERMAL  DIFF  HEAT 

COND’ 

WRITEC*.*)’  NO.  CM  SPACE-CM  CM**2/MIN' 

&  , ’  CAL/MIN-CM-K’ 

DO  99956  J4-l,NOMATL 

WRITEC*. 210) J4, THKCJ4) .SFRQCJ4) . 

&  ALPHCJ4),FKCJ4) 

99956  CONTINUE 
WRITEC*.*)’  ’ 

WRITEC*,*)’  INPUT  BOTTOM  BOUNDRY  DATA' 

WRITEC*,*)'  ' 

IFC.NOT. CLFLUXY.EQ.O))  GO  TO  99958 
WRITEC*, 90)LFLUXY 
WRITEC*, 92)DPRM0- 273. 15 
WRITEC*,*)'  ' 

WRITEC*,*)'  ’ 

GO  TO  99959 

99958  IFC .NOT. CLFLUXY.lt. 0))  GO  TO  99957 
WRITEC*. 90)LFLUXY 

WRITEC*, 95)DPRM1*697. 6 
WRITEC*,*)'  ' 

GO  TO  99959 

99957  WRITEC*, 90) LFLUXY 
WRITEC*, 95)DPRM1*697. 6 

WRITEC*,*)'  ---BOTTOM  SURFACE .  SURFACE  BENEATH  AIRSPACE 

TEMP' 

WRITEC*,*)'  EMISSIVITY  GEOM  SHAPE  EMISSIVITY  GEOM  SHAPE  DEG 
C' 

WRITEC*,*)'  FACTCO.-l.)  FACTCO.-l.)’ 

WRITEC*, 97)BEP, BK, REP. RK,TR-273. 15 
WRITEC*,*)'  ' 

99959  WRITEC*,*)'  ' 

IFCIVEG.EQ.O)  GO  TO  199988 
WRITEC*.*)'  VEGETATION  PARAMETERS' 

WRITEC*,*)'  ' 

WRITEC*,*)'  COVERAGE  STATE  EMISSIVITY  ABSORBTIVITY 
5.  FOLIAGE  HEIGHT' 

WRITEC*,*)  (.0.0  -1.0)  CO.O  -1.0)  CO. 0-1.0) 

«.  CCM) ' 

WRITEC*. 240)SIGF, STATE. EPF, FOIA.HFOL 
WRITEC*,*)'  ' 

WRITEC*,*)'  ' 

WRITEC*.*)'  ' 

GO  TO  199988 

C . 

99985  CONTINUE 

C  TO  CALCULATE-INSOLATION-ON-SLOPE-SURFACE 

C 

C 

C  SOLVE- SOLAR- ZENITH 

ASSIGN  99951  TO  199952 
GO  TO  99952 


c 

C  SOLVE- SOLAR- AZIMUTH 

99951  ASSIGN  99949  TO  199950 
GO  TO  99950 

C 

C  CALCULATE -  SLOPE - ATMOS  - ATTEN- AND - CLOUD - ADJUSTMENTS 

99949  ASSIGN  99947  TO  199948 
GO  TO  99948 
C 

99947  CONTINUE 
C 

99953  GO  TO  199985 

. . 

99955  CONTINUE 
C  TO  ZERO -VARIABLES 
C 

I-O 

GO  TO  199955 

. . 

99952  CONTINUE 

C  TO  SOLVE -SOLAR -ZENITH 
C 

TYME-AMOD ( TIME , 24 . 0 ) 

T0-2.0*PI*(DAY-I.0)/365.0 

DECL-O. 006918 -0.399912*COS(TO)+O.O70257*SIN(T0) 

&  - 0 . 0067  58*COS ( 2 . 0*T0 ) +0 . 000907*S INC  2 . 0*T0 ) 

&  -  0 . 002697*COS ( 3 . 0*T0)+0 . 001480*SIN ( 3 . 0*T0 ) 

ELF-(LAT/180*PI) 

TIMER-(TYME/12*PI)+PI 

IFCTIMER . GT . 2 . *PI)TIMER-TIMER- 2 . *PI 

AA-COS ( DECL) *C0S ( ELF)*COS (TIMER) 

BB-SIN(DECL)*SIN(ELF) 

C-AA+BB 

Z-ACOS(C) 

GO  TO  199952 

. . 


99950  CONTINUE 

C  TO  SOLVE- SOLAR-AZIMUTH 

C 

XNUM—  COS  ( DECL)  *S  IN  (TIMER) 

XDNOM-COS (ELF)*SIN(DECL)-SIN(ELF)*COS(TIMER) 

SAZ-ATAN (XNUM/XDNOM) 

IF(.NOT. (XNUM. LT. 0.0. AND. XDNOM.GT. 0.0))  GO  TO  99944 

SAZ-SAZ+PI 

GO  TO  99945 

99944  IF( .NOT. (XNUM. GT. 0.0. AND. XDNOM.GT. 0.0))  GO  TO  99943 
SAZ-SAZ-PI 

99943  CONTINUE 

99945  GO  TO  199950 

- - - - - - - - - 

99948  CONTINUE 

C  TO  CALCULATE-SLOPE-ATMOS-ATTEM-AND-CLOUD-ADJUSTMENTS 
C 

SICF-COS(Z)*COS(SLOPE)+SIN(Z)*SIN(SLOPE) 

&  *COS(SAZ-SURFAC) 

IF( .NOT. (SICF.LT.O.O.OR.COS(Z) .LE.0.0))  GO  TO  99941 
SUN-0 . 0 
GO  TO  99942 

99941  M-1/C0S(Z) 

IF( .NOT. (M.GE.0.0))  GO  TO  99939 
TAI^O.  02023 

IF(DAY.GE.92.0  .AND.  DAY . LE. 152 .0)TAL--0 . 02290 
TD-5352 . 2/(21 .4-AL0G(RH*ESAT(TA) ) ) 
WATER-EXP(0.07074*(TD-273.15)+TAL) 

AB-0 . 2  7 1*  (  WATER*M  )  **0 . 30  3 

AO-0 . 085  -  0 . 247*AL0G10 ( PRESS/1000 . *1 . /M) 

ARG1-((1. -AB)*0.349+(1. -A0)/(1.-a0*0.2)*0.651) 

GO  TO  99940 

99939  ARGl-1.0 

99940  qP-2 . 0*ARG1 
QO-QP*SICF 

IF(.NOT. (NCLOUD.EQ.O))  GO  TO  99937 

SUN-QO 

GO  TO  99938 

99937  CONTINUE 
ARG2--(BCL(NCLOUD)- .059)*M 
CTF-(ACL(NCLOUD)/94.4)*EXP(ARG2) 

SUN-qO- ( (CLOUD*CLOUD)*(qO-qO*CTF) ) 

99938  CONTINUE 

99942  CONTINUE 

GO  TO  199948 

C . 

99983  CONTINUE 

C  TO  CALCULATE -TABLE -SLOPE -AND -INTERCEPT 
I-I 

GO  TO  99935 

99936  IF(I.GT.IO)  GO  TO  99934 
99935  IMAX-MAX(I) 

IF(IMAX.Eq.O)GO  TO  99931 
J-1 

GO  TO  99932 

99933  IF(J.Eq.IMAX)  GO  TO  99931 

99932  FMM(J , I)-(YYY(J+1 ,1) -YYY(J . I) )/(XXX(J+l , I) -XXX(J , I) ) 
BBB(J , I)-YYY(J , I) -FMM(J ,I)*XXX(J .1) 


O  Ci 


J-J+1 

GO  TO 

99933 

99931 

I-I+l 

GO  TO 

99936 

99934 

GO  TO 

199983 

c - *  - 

99930 

CONTINUE 

TO  GET-TABLE-VALUES 

IMAX-MAX(NTABL) 

IJ-1 

IF( .NOT. (XN.GE.XXX(IMAX,NTABL) ))  GO  TO  99928 
YN-YYY ( IMAX , NTABL) 

GO  TO  99929 

99928  IF(IJ.EQ.IMAX+1)  GO  TO  99927 
JJ-IJ 

IF(. NOT. (XXX(IJ, NTABL). LT.XN))  GO  TO  99923 

IJ-IJ+1 

GO  TO  99928 

99925  IF( .NOT. (XXX(IJ, NTABL) .EQ.XN))  GO  TO  99924 
YN-YYY(JJ, NTABL) 

IJ-IMAX+1 
GO  TO  99928 

99924  IF(. NOT. (XXX(IJ. NTABL). GT.XN))  GO  TO  99923 
JJ-JJ-1 

YN-FMMUJ  ,NTABL)*XN+BBB(JJ  .NTABL) 

IJ-IMAX+1 
99923  CONTINUE 

GO  TO  99928 
99927  CONTINUE 

99929  CO  TO  199930 

. . 

99981  CONTINUE 

C  TO  CALCULATE-SURFACE-AND-LAYER-TEMPERATURE 

C 

C  SET-UP-INITIAL-CONDITIONS 

assign  79921  TO  199922 
GO  TO  94922 
C 

C  PRINT-OUTPUT-HEADING 

99921  ASSIGN  99919  TO  199920 
GO  TO  99920 
C 

99919  CONTINUE 
C 

C  RUN -HEAT -FLOW -PROGRAM 

ASSIGN  99917  TO  199918 
GO  TO  99918 
C 


C  SET-UP-AND-PRINT-OUTPUT 

99917  ASSIGN  99915  TO  199916 
GO  TO  99916 
C 

99915  IFdRETRN.EQ.  1)  GO  TO  4100 
GO  TO  99919 
4100  CONTINUE 

GO  TO  199981 

C- . 

99922  CONTINUE 

C  TO  SET-UP- INITIAL-CONDITIONS 
C 

PTVME-TOTTIM-24,0 

TIME-XXX(l.l) 

DIST-0. 

I FLAG-0 

IF  (TFRQ.LE.O)  TFRQ-TOTTIM 
DELT-TFRQ/60. 

ITIME-MAXl (TOTTIM/DELT+ .9,1.1) 

IX-1 

IY-1 

GO  TO  99913 

99914  IF(IX.GT.NOMATL)  GO  TO  99912 
99913  INTR(IX)-IY 

IF  (SFRQ(IX) .LE.O. )  SFRQ(IX)-THK(IX)/10 . 

NX  ( IX) -MAXI  ( THK  ( IX ) /SFRQ  ( IX)  + .  9 , 1 . 1 ) 
c  wrice(9 ,*) ' sfrq(ix) .nx(ix)-  ' , sfrq( ix) , nx( ix) 

RR(IX)-60.0*DELT/(SFRQ(IX)*SFRQ(IX)) 
stabn(ix)-alph(ix)*rr( lx) 
c  uri te (9 , *) ' stabn( ix)-  '.sCabn(ix) 

if (stabn( ix) . gc . 0 . 5) then 
sCabn(ix)-0. 5 

c  write (9 , *) ' stabn( ix)-  ',stabn(ix) 

end  if 
LAYERS-0 
GO  TO  99910 

99911  IF(LAYERS.GT.NX(IX))  CO  TO  99909 
99910  XN-DIST 
NTABL-5 
DEPTH(IY)-XN 
C 

C  GET-TABLE-VALUES 

ASSIGN  99908  TO  199930 
GO  TO  99930 
C 

99908  TEMP-YN 

STOR(l,IY)-TEMP 

STOR(5,IY)-TEMP 

STOR(6,IY)-FK(IX) 

STOR(7.IY)-RHOC(IX) 

ST0R(4,IY)-0. 

STOR(2,IY)-FK(IX) 

STOR(3,IY)-RNOC(IX) 

c  write(9 , *) ' STOR  MATRIX-  ' , stor( 1 . iy) . scor ( 5 , iy) . stor (6  ,  iy) , 

c  &  stor ( 7 , iy) , stor (4 , iy) , stor(2 , iy) , stor( 3 , iy) 


IY-IY+1 

DIST-DIST+SFRQ(IX) 

LAYERS-LAYERS+1 
GO  TO  99911 
99909  IX-IX+1 

DIST-DIST-SFRQ(IX-l) 

GO  TO  99914 
99912  JMAX-IY-1 

1NTR(IX)-JMAX 

NPRNT-MAXl ( TPRNT/TFRQ+ .9,1.1) 

IPRNT-NPRNT 
GO  TO  199922 

C . - . 

99920  CONTINUE 

C  TO  PRINT-OUTPUT- HEADING 

C 

IF(OUTCD.EQ.l)GO  TO  1610 
IF(IVEG.GT.O)  GO  TO  1420 
WRITE(*.350) 

WRITE(*,360) 

URITE(*,370) 

WRITE(*,380) 

GO  TO  199920 
1420  WRITE(*,310) 

WRITE(*,320) 

WRITE(*,330) 

WRITE(*, 340) 

GO  TO  199920 
1610  WRITE(*,390) 

GO  TO  199920 

C . 

99918  CONTINUE 

C  TO  RUN -HEAT -FLOW -PROGRAM 

C 

c  write (9 ,*) ' lflag2-  '.iflag2 

IF(.NOT.(IFLAG2.EQ.O))  GO  TO  99907 
ITER-NIT 

c  write(9 ,*)' iter , time , delt-  iter , time  ,  del t 

TIME-TIME+DELT 
99907  ZZA-STOR(5,l) 

ZZB-STOR(5.JMAX) 

TEML-ZZA 

TEMR-ZZB 

C 

C  CALCULATE-BOUNDRY-CONDITIONS  calculate  energy  budget  te 

IF(IVEG.EQ.0)GO  TO  930 
ASSIGN  99905  TO  199800 
GO  TO  99800 

930  ASSIGN  99905  TO  199906 
GO  TO  99906 
C 


C  CALCULATE-UPPER-BOUNDARY-VALUES  using  energy  budget,  find  top 

tmp 

99905  IF(IVEG.EQ.O)  GO  TO  900 
ASSIGN  99903  TO  199797 
GO  TO  99797 

900  ASSIGN  99903  TO  199904 
GO  TO  99904 
C 

99903  IX-1 
J-2 

C  ***  NOMAT  CHANGED  TO  NOMATL  NOV  22  85  *** 

IMATL-NOMATL 

IF(. NOT. (NOMATL. EQ.l))  GO  TO  99901 
I2-NX(IX)-1 

IF( .NOT. (IZ.GT.O))  GO  TO  99900 
C 

C  CALCULATE-INSIDE-MATERIAL-VALUES 

ASSIGN  99898  TO  199899 
GO  TO  99899 
C 

99898  CONTINUE 

99900  GO  TO  99902 

99901  GO  TO  99896 

99897  IF(IMATL.LT. 1)  GO  TO  99895 
99896  IF( .NOT. (IMATL.GT.l))  GO  TO  99893 
IZ-NX(IX)-1 

c  wrice(9,*)'**FLAG  51**iz-Mz 

IF( .NOT. (IZ.GT.O))  GO  TO  99892 
C 

C  CALCULATE- INSIDE  MATERIAL-VALUES 

ASSIGN  99891  TO  199899 
GO  TO  99899 
C 

99891  CONTINUE 
C 

C  CALCULATE- INTERFACE- VALUES 

99892  ASSIGN  99889  TO  199890 
GO  TO  99890 

99889  GO  TO  99894 

99893  IZ-NX(IX)-1 

IF( .NOT. (IZ.GT.O))  GO  TO  99888 
C 

C  CALCULATE- INSIDE-MATERIAL-VALUES 

ASSIGN  99887  TO  199899 
C 

GO  TO  99899 

99887  CONTINUE 

99888  CONTINUE 

99894  IMATI^IMATL-1 
GO  TO  99897 

99895  CONTINUE 
C 

C  CALCJLATE -  LOWER -  BOUNDARY - VALUES 

99902  ASSIGN  99883  TO  199886 
GO  TO  99886 

C 

99883  GO  TO  199918 

C . 


99906  CONTINUE 

C  TO  CALCULATE-BOUNDRY-CONDITIONS 
B  -  -FK(1) 

T-TIME 

C 

C  CALCULATE -  BOTTOM - BOUNDRY - HEAT - TERMS -APRM-DPRM-BPRM 

ASSIGN  99880  TO  199881 
GO  TO  99881 
C 
C 

C  ATMOSPHERIC -INFRARED -EMISSION-ATERM 

99880  ASSIGN  99878  TO  199879 
GO  TO  99879 

C 

99878  CONTINUE 

C  CALCULATE- SOLAR- INSOLATION- FOR-DAY-AND-TIME 
DAY-XXX(1,8) 

ASSIGN  42  TO  199985 
GO  TO  99985 
42  BTERM-SUN 

I F ( BTERM . GT . 0 ) BTERM-BTERM*SMALLA 
C 

C  CALCULATE- CONVECTION -HTERM 

ASSIGN  99876  TO  199877 
GO  TO  99877 
C 
C 

C  CALCULATE  -  EVAPORATIVE  - HEAT  -  LOSS  - DTERM 

99876  ASSIGN  99874  TO  199875 

GO  TO  99875 
C 

99874  D  -  ATERM  +  BTERM  -  HTERM-DTERM 
521  CONTINUE 

GO  TO  199906 

C . 

99881  CONTINUE 

C  TO  CALCULATE-BOTTOM-BOUNDRY-HEAT-TERMS-APRM-DPRM-BPRM 

c 

BPRM-TB 

IF( .NOT. (TB.EQ.0.0))  GO  TO  99872 

APRM-1.0 

DPRM-DPRMO 

GO  TO  99873 

99872  APRM-FACTE*TEMR*TEMR*TEMR 
DPRM-DPRMl 

99873  GO  TO  199881 

C . 

99879  CONTINUE 

C  TO  ATMOSPHERIC-INFRARED-EMISSION-ATERM 
XN-TIME 
NTABL-2 
C 


C  GET-TABLE-VALUES 

ASSIGN  99867  TO  199930 
GO  TO  99930 
C 

99867  RH-YN 

XN-TIME 

NTABI^IO 

C 

C  GET -TABLE -VALUES 

ASSIGN  99866  TO  199930 
GO  TO  99930 
C 

99866  TA-YN 

XN-TIME 

NTABL-3 

C 

C  GET-TABLE-VALUES 

ASSIGN  99865  TO  199930 
GO  TO  99930 
C 

99865  CLOUD- YN 
TAK-TA 

TAC-(TAK-273.15) 

EA-6 . 108*RH*EXP( (AC*TAC)/(TAK-BC) ) 

ALPHI-(0.61+0.05*SQRT(EA))*(l,0+(CLR(NCLOUD)*fCLOUD**2) ) ) 
DOWNIR-0 . 8132E- 10*TAK**4*ALPH1 
ATERM-DOWNIR 
GO  TO  199879 

C . 

99877  CONTINUE 

C  TO  CALCULATE-CONVECTION-HTERM 
C 

XN-TIME 

NTABL-6 

C 

C  GET -TABLE -VALUES 

ASSIGN  99862  TO  199930 
GO  TO  99930 
C 

99862  SPEED-YN 
TAK-TA 
ZASH-ZA 
TSK-TEML 

RHOA--0 . 001*0 . 348*PRESS/TAK 
1200  THETAZ-TAK*FACTH 
THETAS-TSK*FACTH 
DTHETA- ( THETAZ -  THETAS ) /ZASH 
DU-SPEED/ZASH 

THETAV-(THETAZ+THETAS)/2 . 0 
RI-G*DTHETA/ (THETAV*DU**2 ) 

COEl-15.0 

COE2-1.175 

EX-.75 


IF(TSK.GT,TAK)GO  TO  31 
IF(RI.GT.n.2)RI-. 19999 
COEl-5 . 0 
COE2-1 . 0 
EX-2.0 

31  HTER-RHOA*KSQ*2ASH**2*DU 
&  *(COE2*(1.0-COE1*RI)**EX) 

HTERM-HTER*CP*DTHETA 

c  write(9 , *) ' tak, zash, Csk, rhoa-  ' , Cak.zash, tsk, rhoa 

c  wrlte(9 ,*)' thetaz , thetas .dtheta-  thetaz , thetas , dtheta 

c  write(9 ,*) 'du, thetav, ri-  ' ,du, thetav, ri 

c  write(9 ,*) 'hter .hterm-  ' ,hter,hterm 

c  write(9 ,*) 'ksq , coe2 ,coel ,ex-  ' ,ksq , coe2 , coel , ex 

99864  GO  TO  199877 

C . 

99875  CONTINUE 

C  TO  CALCULATE-EVAPORATIVE-HEAT-LOSS-DTERM 
C 

IF( .NOT. (TEML.GT.TA))  GO  TO  99860 
XN-TIME 
NTABL-2 
C 

C  GET-TABLE-VALUES 

ASSIGN  99859  TO  199930 
GO  TO  99930 
C 

99859  RH-YN 
XN-TIME 
NTABL-10 

C 

C  GET -TABLE -VALUES 

ASSIGN  99858  TO  199930 
GO  TO  99930 
C 

99858  CTEMA-YN 

KTEMPA-CTEMA 

CTEMA-CTEMA-273.15 

KTEMPG-TEML 

ES-EXP((AC*(KTEMPG-273. 15))/(KTEMPC-BC))*6, 1071 

EA-EXP( (AC*CTEMA)/(KTEMPA-BC) )*6 . 1071*RH 

DG-0 . 622/PRESS*(EA- ES)*WET/ZA 

XL-597 . 3-0 . 566*(CTEMA+KTEMPG-273 . 15)/2 .0 

DTERM-HTER*XL*DG 

GO  TO  99861 

99860  DTERM-0.0 

99861  GO  TO  199875 

C . . 


99904  CONTINUE 

C  TO  CALCULATE-UPPER-BOUNDARY-VALOES 

Tl-stabn(X)*(STOR(l,3)-2.*STOR(1.2)+STOR(l,l;>+STOR(l,2) 
c  write(9 , *) ' storl , 2 , 3-  ' , scor(l , 1) , stor(l , 2) , scor (1 , 3) 

c  write(9,*) 'alph(l) ,rr(l)-  ' , alph(l) , rr (1) 

III-O 

830  III-III+l 

c  write(9 , *) ’ iii-  '.iii 

T2-STOR ( 5 , 1) **4*FACTA*S  FRQ ( 1 ) +FK ( 1 ) *STOR ( 5 , 1 ) 

&  -(FK(1)*T1+D*SFRQ(1)) 

c  write(9 ,*)' facta, sfrq(l) , fk(l)-  facta , sftq( i) , fk( 1) 

c  wrlte(9 ,*) ' tl , d, stor3a-  ' , tl ,d. sCor(5 , 1) 

T2-T2/(4.*FACTA*SFRQ(l)*STOR(5.1)**3+FK(l) -SFRQ(1)*DDDT) 
3TOR(5,l)-STOR(5,l)-T2 
c  write(9 ,*) ' stor4 , t2-  ' , stor (5 , 1) , t2 

TEML-ST0R(5,1) 

GTERM--FK(1)*(ST0R(5,1)-T1)/SFRQ(1) 

ASSIGN  825  TO  199877 
GO  TO  99877 

825  ASSIGN  810  TO  199875 
GO  TO  99875 

810  DNEW-ATERM+BTERM-HTERK-DTERM 

IF(ABS(T2) .LT. 0.005  .OR.  Ill .GT. 30)GO  TO  199904 
DDDT— (DNEW-D)/T2 
D-DNEW 
GO  TO  830 
C 

C . 

99899  CONTINUE 

C  TO  CALCULATE- INSIDE-MATERIAL-VALUES 
C 

GO  TO  99856 

99857  IF(IZ.LE.O)  GO  TO  99855 
99856  CONTINUE 

STOR(5,J)-STOR(l,J)+stabn(lx)*(STOR(l,J-l)-2.*STOR(l ,J) 

&  -tSTORd.J+D) 

J-0>1 
IZ-IZ-1 
GO  TO  99857 
99855  GO  TC  199899 

C . 

C  TO  CALCULATE -INTERFACE -VALUES 
C 

99890  CONTINUE 

BCOEF-STOR ( 6 , J - 1 )/S FRQ( IX ) 

DCOEF-STOR( 6 , J+1 ) /SFRQ( IX+1 ) 
c  write(9,*) 'stor9, 10-  ' , scor(6 , ) - 1) . stor(6 , j  +  1 ) 

CCOEF-BCOEF+DCOEF 

ACOEF-BCOEF/(2.*stabn(ix))+DCOEF/(2.*stabn(ix+l) ) 


STOR( 5 . J ) -STOR ( 1 , J ) + (BCOEF*STOR( 1 , J - 1 ) - CCOEF*STOR (1 , J ) +DCOEF* 
&  ST0R(1,J+2))/AC0EF 

c  write(9 ,*) ’ storll , 12 , 13-  ' , stor(l , j ) , stor (1 , j - 1) , scor (1 , j +2) 

c  write(9 ,*) ' storl4-  ’,stor(5,J) 

STOR( 5 , J+1) -STOR ( 5 , J ) 

IX-IX+1 

J-J+2 

GO  TO  199890 

C . 

C  TO  CALCULATE -LOWER -BOUNDARY -VALUES 
99886  IF(LFLUXY.EQ.O)  GO  TO  ?80 
I-l 

R2-FACTD 
870  CONTINUE 

Rl-S IGMA*BEP*BK*STOR ( 5 , J ) **4 

Gl— FK(NOMATL)*(STOR(5,J)-STOR(l,J-l))/SFRQ(NOMATL) 
F2-4.0*SIGMA*BEP*BK*STOR(5,J)**3-FK(NOMATL)/SFRQ(NOMATL) 

F2-  -  (R2-R1+G1+DPRM)/F2 
STOR(5,J)-STOR(5,J)  +  F2 
c  wrlce(9 ,*) ' storl5-  ',stor(5,j) 

I-I+l 

IF(ABS(F2) .GT.0.01  .AND,  I.LE.30)  GO  TO  870 
880  IF(LFLUXY.EQ.O)  STOR( 5 , J)-ST0R(5 , J ) 

GO  TO  199886 

C . 

99916  CONTINUE 

C  TO  SET-UP-AND-PRINT-OUTPUT 
C 

IFLAG2-0 

IRETRN-0 

IF(ITER.LE.l)  GO  TO  1243 
ITER-ITER- 1 

IF  (lEFSWT.NE.O)  GO  TO  1328 
IF  (IPRNT.LE.l)  GO  TO  1244 
IF  (ITIME.GT.l)  GO  TO  1328 

1244  CONTINUE 
C 

C  PRINT -OUTPUT 

ASSIGN  99846  TO  199847 
GO  TO  99847 
C 

99846  GO  TO  1328 

1245  IF  (ITIME.GT.l)  GO  TO  1269 
C 

C  PRINT-OUTPUT 

ASSIGN  99845  TO  199847 
GO  TO  99347 
C 


99845  IRETRN-1 

GO  TO  1335 
1269  ITIME-lTIME-1 

IF  (IPRNT.LE.l)  GO  TO  1279 
IPRNT-IPRNT-1 
GO  TO  1303 
1279  CONTINUE 
C 

C  PRINT -OUTPUT 

ASSIGN  99844  TO  199847 
GO  TO  99847 
C 

99844  IPRNT-NPRNT 
1303  J-1 

IZ-NOMATL 
1306  IX-NX(IZ)+1 
1311  STOR(l,J)-STOR(5,J) 

STOR(2.J)-STOR(6,J) 

STOR(3.J)-STOR(7,J) 

J*J+1 

IF  (IX.LE.l)  GO  TO  1329 

IX-IX-1 

GO  TO  1311 

1329  IF(IZ.LE.l)  GO  TO  1335 
IZ  -  lZ-1 
GO  TO  1306 
1328  IF1AG2-1 
1335  CONTINUE 

GO  TO  199916 

C . 

99847  CONTINUE 
C  TO  PRINT-OUTPUT 
C 

TFACK-1 . 0 

IF(. NOT. (TIME. GT.PTYME))  GO  TO  99843 
DO  99842  JKX-l.NOMATL-t-l 
IJ-INTR(JKX) 

TITLE(JKX)-(STOR(5,IJ)-273.15) 

99842  CONTINUE 
NDX-TIME 

IF(NDX.Eq.0)NDX-l 
bottp-stor(5 , j ) -273 . 15 
klO-klO+1 

IF(IVEG.Eq.l)  GO  TO  1110 
IGBR-5.67E-8*EPSN*STOR(5,l)**4 
c  wrlCe(9 ,*)'**- IGBR-epsn, stor (5 , 1)-  ' , epsn . stor ( 5 , 1) 

ISOL-BTERM/SMALLA*697 .6+0.5 
lABSOR-I  SOL*SMALU 
IATERM-ATERM+697 . 6 
IHTERM-HTERM+697 .6+0.5 
IDTERM-DTERM+697 .6+0.5 
1110  continue 


c  write (9 ,*)' Bottom  temp.-  ',stor(5.j) 

datmx(klO , l)-time 
datmxCklO ,2)-title(l) 

ASSIGN  1400  TO  11410 
if(iveg.eq. 1)  GO  TO  1410 
1400  continue 

if(iveg.eq.O)  then 

WRITE(9 , 190)AMOD(TIME, 24 . ) .TITLE(l) , IGBR, ISOL, lABSOR 
&  lATERM , IHTERM , IDTERM 

datmx(klO , 3)-igbr 
datna(klO  ,4)-isol 
datmx(klO , 5)-iabsor 
datmx(klO , 6)-iaterm 
datmx(klO, 7)-ihterm 
datmx(klO , 8)-idterm 
datmx(klO , 9)-bottp 

else 

WRIT£(9,270)AM0D(TIME.24.) .ISURFG+IREFRA.ISURFG.TEFFR 

273.15, 

&  TEFF-273.15,TEML-273.15.TF-273.15, ISOL 

datmx(klO , 3)-rlu*697 . 6 
datmx(klO ,4)-sol*697 . 6 
datmx(klO , 5)-smalla*sg*697 . 6 
datmx(klO , 6)-rld*697 . 6 
datmx(klO , 7)-hsg*697 . 6 
datmx(klO, 8)-elg*xll*697 . 6 
datmx(klO , 9)-bottp 
endif 

IF(AN.EQ. 'Y' )THEN 

WRITE<6 , 190)TIME.TITLE(1) , IGBR. ISOL. lABSOR, 

&  lATERM, IHTERM, I DTERU 
c  write(9 ,*)' Bottom  temp.-  '.stor(5,j) 

ENDIF 

PRINT  *,'TIME,  AIR  TEMP  '.TIME,  TA-273.I5 
99843  GO  TO  199847 

C . 

99800  CONTINUE 

C  TO  CALCULATE-BOUNDARY-CONDITIONS-WITH-VEG 

T-TIME 
XN-TIME 
NTABL-6 
C 

C  GET-TABEL-VALUES 

ASSIGN  960  TO  199930 
GO  TO  99930 
C 

960  UA-YN 

XN-TIME 

NTABI^4 

C 

C  ATMOSPHERIC- INFRARED- EMISSION-ATERM 

ASSIGN  980  TO  199879 
GO  TO  99879 
C 


C  Calculate  solar 
980  CONTINUE 

DAY-XXX(1.8) 
assign  43  TO  199985 
GO  TO  99985 
43  SOL-SUN 

C 

IF(UA.LT.10.0)UA-10.0 

UAF-O . 83*SIGF*UA*SQRT(CHH)+(1 . -SIGF)*UA 
DELTMP-5 . 

CF-0.01*(1.+30.0/UAF) 

DU-(UA-UAF)/ZA 
RS-1/(.05+.0021*(SOL*697  t.)) 

RC-RS*STATE/(7.0*SIGF) 

ATF(1)-TF 

ASSIGN  1210  TO  1950 
GO  TO  950 
1210  CONTINUE 

FEB(1)-FENB 

NDEX-0 

1240  TF-TF+DELTMP 

NDEX-NDEX+1 
ASSIGN  1220  TO  1950 
GO  TO  950 
1220  CONTINUE 

FEB(2)-FENB 

IF(FEB<1)*FEB<2) .LT.0,0)  GO  TO  1230 
IF<ABS(FEB(2)) .GT.A3S(FEB(1)))DELTMP— 5. 

IF(NDEX,LT.100)GO  TO  1240 

WRITE(*,*) 'FOLIAGE  ENERGY  BUDGET  HAS  NOT  CROSSED  X-AXIS' 
WRITE (*,*) 'AFTER  100  SEARCH  STEPS.  CHECK  INPUT  DATA.' 
STOP 

1230  CONTINUE 

ATF(2)-TF 

1270  SL0PE1-(FEB(2)-FEB(1))/(ATF(2)-ATF(1)) 

BINT-FEB(l) -SL0PE1*ATF(1) 

TF0--BINT/SL0PE1 

IF(ABS(TF-TF0) .LE.0.001)GO  TO  1260 
TF-TFO 

ASSIGN  1250  TO  1950 
GO  TO  950 
1250  CONTINUE 

IF(FENB*FEB(2) .GT.C.0)lP-2 

IF(FENB*FEB(1) .GT.0.0)IP-1 

ATF(IP)-TF 

FEB(IP)-FENB 

GO  TO  1270 

1260  GO  TO  199800 

. . 


C  TO  CALCULATE-UPPER-BOUNDARY-VALUES-FOR-FOLAGE 

99797  CONTINUE 

DELTMP-5 . 

ATF(1)-TEML 
ASSIGN  1310  TO  11300 
GO  TO  1300 
1310  CONTINUE 

FEB(1)-FENB 

NDEX-0 

1340  TEML-TEML+DELTMP 

NDEX-NDEX+1 
ASSIGN  1320  TO  11300 
GO  TO  1300 
1320  CONTINUE 

FEB(2)-FENB 

IF(FEB(1)*FEB(2).LT.0.0)  GO  TO  1330 
IF(ABS(FEB(2)) .GT.ABS(FEB(l)))DELTMP--5 . 
IF(NDEX.LT.100)GO  TO  1340 

WRITE (*,*) 'GROUND  ENERGY  BUDGET  HAS  NOT  CROSSED  X-AXIS' 
WRITE(*,*) 'AFTER  100  SEARCH  STEPS.  CHECK  INPUT  DATA  ' 
STOP 

1330  CONTINUE 

ATF(2)-TEML 

1370  SL0PE1-(FEB(2)-FEB(1))/(ATF(2)-ATF(1)) 

BINT-FEB(1)-SLX)PE1*ATF(1) 

TFO--BINT/SLOPE1 

IF(ABS(TEML-TFO)  .LE.O.OODGO  TO  1360 
TEML-TFO 

ASSIGN  1350  TO  11300 
GO  TO  1300 
1350  CONTINUE 

IF(FENB*FEB(2) ,GT.0.0)IP-2 
1F(FENB*FEB(1) ,GT.0.0)IP-1 
ATF(IP)-TEML 
FEB(IP)-FENB 
GO  TO  1370 

1360  STOR(5.1)-TEML 

GO  TO  199797 

C . 

C  TO  CALCULATE- ENERGY- BUDGET 
950  TAF-(1.  -SIGF)*TA-^SIGF*(0.3*TA■^0.6*TF■^0,1*TEML) 

DTHETA- ( TA - TF ) *FACTH/ZA 
THETA V-  ( TA-t-TF)  *FACTH/2 . 0 
RI-G*DTHETA/ (THETA V*DU**2) 

RHOAF--0.001*. 348*PRESS/((TF+TA)/2. ) 

COEl-15. 

COE2-1.175 

EX-.75 

IF(RI.LE.O.)GO  TO  1280 

IF(RI.GT.0.2)RI-0.199 

COEl-5 

COE2-1. 

EX-2.0 


1280  CONTINUE 

HTER-RH0AF*KSQ*ZA*-*2*DU 
&  *C0E2*(1. -C0E1*RI)**EX 

C  HSF-1.1*7.*SIGF*CP*CF*UaF*(TF-TAF)*60. 

HSF-HTER*CP*DTHETA*60 . 

XL-597.3-0.566*TAF 

RA- ( ALOG ( ( ZA- ZDSP ) /ZO ) *COE2* ( ( 1 . - C0E1*RI ) **E<) ) ** 2 
&  /(.16*UA) 

RDP-RA/(RS+RA) 

QF-RDP*QSAT(TF)+(1 . -RDP)*QAF 

QAF-(1 .  -SIGF)*Q(TA)+SIGF*(Q(TA)*0 . 3-K!F*0 . 6+00*0 . 1 ) 

EF— (RHOAF+CP/O . 66)*(ESAT(TF) ’E(TA) )/(RA+RC)*60. 

IF(EF.LT.0.0)EF-0.0 

SHRW-F01-A*SOL 

XLNGW-EPF*ATERM 

TGA-EPF*EPSN/EPl*SIGMA*TEiML**4 

TF4-(EP1+EPSN)/EP1*EPF*SIGMA*TF**4 

FENB-S IGF* ( SHRW+XLNGU+TG4 - TF4 ) - HS  F - EF 

GO  TO  1950 

C . 

C  TO  CALCULATE-ENERGY-BUDGET-FOR-GROUND 

1300  CONTINUE 

T1-ALPH(1)*RR(1)*(ST0R(1,3) -2.*STOR(l,2)+STOR(l.l) ) 

&  +ST0R(1,2) 

TF4-S1GMA*TF**4 

TG4-SIGMA*TEML**4 

QG-WET*QSAT(TEML)+(1. -WET)*QAF 

RHOAG-0 . 001*0 . 348*PRESS/TAF 

XLl-597. 3-0. 565*(TAF+TEML-2. 0*273. 15)/2. 

SG-(1.-S1GF)*S0L 

RLU-(1, •SIGF)*(EPSN*TG4+(l. -EPSN)*ATERM) 

&  +SIGF*(EPSN*TG4+(1. -EPSN)*EPF*TF4)/EP1 

RLD-d.  ■SIGF)*ATERH+SIGF*(EPF*TF4+(1.  -EPF)*EPSN*TG4)/E?1 
HSG-RHOAC*CP*CHG*UAF*(TEMI,-TAF)*60. 
ELG-RHOAG*CHG*UAF*(QG-QAF)*60 . 

FENB-SMALLA*SG  -  RLU+RLD  -  HSG  -  ELG*XU  *  ( T1  -  TEML) /S  FRQ  ( 1 )  *FK  ( 1 
GO  TO  11300 

. . 

C  TO  CAiCUTATE-RADlANCE-VALUES 
1410  CONTINUE 

REFRAD-((1. -SIGF)*(l-EPSN)+SIGF*(l-EPF))*DOWNIR*697 . 6 

FOLGB-EPF*5 . 67E- 8*TF**4 

GRNt)GB-EPSN*5 .67E-8*TEMI.**4 

SURFGB-SIGF*FOLGB+(l, -SIGF)*GRNDGB 

EEF-SIGF*EPF+(1. -SIGF)*EPSN 

TEFF- ( SimFGB/ 5 , 6 7E  -  8  )  **  .  2 5 

ISURFG-SURFGB+ . 5 

TEFFR- ( ( SURFGB+REFRAD ) /(5.67E-8))**,25 
IREFRA-REFRAJHO ,  5 
ISOl^S0L*697.6+0.5 
GO  TO  11410 
C 

END 


^  nr 

V',/ 


